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Abstract 

In  this  dissertation,  I  explored  metapopulation  dynamics  and  population  connectivity, 
with  a  focus  on  the  softshell  clam,  Mya  arenaria.  I  first  worked  towards  developing 
a  method  for  using  elemental  signatures  retained  in  the  larval  shell  as  a  tag  of  natal 
habitat.  I  designed  and  implemented  an  experiment  to  determine  whether  existing 
methods  commonly  used  for  fishes  would  be  applicable  to  bivalves.  I  found  that  the 
instrumentation  and  setup  I  used  were  not  able  to  isolate  and  measure  t  he  first  larval 
shell  of  M.  arenaria.  In  concert  with  developing  this  method  for  bivalves,  1  reared 
larval  M.  arenaria  in  the  laboratory  under  controlled  conditions  to  understand  the 
environmental  and  biological  factors  that  may  influence  elemental  signatures  in  shell. 
My  results  show  that,  growth  rate  and  age  have  significant  effects  on  juvenile  shell 
composition,  and  that  temperature  and  salinity  affect  larval  and  juvenile  shell  com¬ 
position  in  variable  ways  depending  on  the  element  evaluated.  I  also  examined  the 
regional  patterns  of  diversity  over  the  current  distribution  of  M.  arenaria  using  the 
mitochondrial  gene,  cytochrome  oxidase  I  (COI).  I  found  minimal  variability  across 
all  populations  sampled,  suggesting  a  recent  population  expansion  in  the  Northwest 
Atlantic.  Finally,  I  employed  theoretical  approaches  to  understand  patch  dynamics  in 
a  two-patch  metapopulation  when  one  patch  is  of  high  quality  and  the  other  low  qual¬ 
ity.  I  developed  a  matrix  metapopulation  model  and  compared  growth  rate  elasticity 
to  patch  parameters  under  variable  migration  scenarios.  I  then  expanded  the  model 
to  include  stochastic  disturbance.  I  found  that  in  many  cases,  the  spatial  distribution 
of  individuals  within  the  metapopulation  affects  whether  growth  rate  is  most  elastic 
to  parameters  in  the  good  or  bad  patch. 
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Chapter  1 


Introduction 

1.1  Metapopulations 

The  concept  of  a  met.apopulation,  first  introduced  by  Levins  (1969;  1970),  has  received 
much  attention  in  marine  ecology  over  the  last  two  decades  (e.g.  Roughgarden  and 
Iwasa,  1986;  Hanski,  1999;  Thorrold  et,  ah,  2001;  Kritzer  and  Sale,  2006).  Original 
definitions  of  the  metapopulation  centered  around  extinction/recolonization  events 
as  the  most  important  process  to  examine  in  metapopulation  studies  (Levins,  1969). 
More  recently,  the  concept  of  a  metapopulation  has  been  expanded  to  include  many 
other  questions  about  interconnected  populations,  including  population  size,  struc¬ 
ture,  and  long-term  trends. 

In  2004,  Kritzer  and  Sale  examined  the  validity  of  Levins’  original  metapopula¬ 
tion  definition  from  the  perspective  of  marine  ecology,  drawing  heavily  from  Hanski 
(1999),  and  re-defined  a  metapopulation  based  on  the  two  spatial  scales  at  which 
metapopulations  exist:  the  patch  (i.e.  subpopulation  or  local  population)  scale,  and 
the  regional  scale  (Kritzer  and  Sale,  2004,  p.  138).  More  recently,  Sale  et  al.  (2006) 
defined  the  degree  of  connection  between  patches  more  explicitly  for  their  definition 
of  a  metapopulation: 

We  define  a  metapopulation  as  a  system  in  which  (1)  local  populations 
inhabit  discrete  habitat  patches  and  (2)  interpatch  dispersal  is  neither 


15 


so  low  as  to  negate  significant  demographic  connectivity  nor  so  high  as 
to  eliminate  any  independence  of  local  population  dynamics,  including  a 
degree  of  asynchrony  with  other  local  populations.  (Sale  et  ah,  2006,  p.  6) 

The  degree  to  which  patches  are  interconnected  is  determined  by  the  migration  of 
individuals  between  those  patches,  which  then  dictates  their  demographic  connectiv¬ 
ity.  Because  of  this  link  between  migration  and  connectivity,  Kritzer  and  Sale  (2004) 
suggest  that  estimates  of  exchange  rates  should  be  a  priority  for  future  research.  Mi¬ 
gration  is,  in  fact,  one  of  the  most  important  features  for  determining  metapopulation 
dynamics  (Stacey  et  al.,  1997).  Its  importance  extends  beyond  the  effects  of  migra¬ 
tion  on  local  population  dynamics.  Migration  facilitates  gene  flow  among  patches 
and  therefore  affects  metapopulation  genetic  structure  (Gillespie,  1981),  and  is  es¬ 
pecially  critical  for  long-term  metapopulation  success  when  patches  are  subjected  to 
stochastic  environmental  variability  (Howe  and  Davis,  1991;  Bascompte  et  ah,  2002; 
Hill  et  ah,  2002). 

Although  the  metapopulation  concept  evolved  primarily  in  the  terrestrial  liter¬ 
ature,  it  is  highly  relevant  to  marine  populations.  This  was  not  recognized  in  the 
literature  any  appreciable  numbers  until  the  mid-90’s  (Grimm  et  ah,  2003).  Prior 
to  that  time,  marine  ecologists  had  focused  on  the  importance  of  larval  recruitment 
for  structuring  marine  populations  (Connell,  1985),  later  recognizing  the  impacts  of 
recruitment  variability  in  population  dynamics  on  coral  reefs  (Sale,  1977)  and  along 
rocky  coastlines  (Caffey,  1985;  Roughgarden  et  ah,  1985).  Although  evidence  from 
larval  dispersal  and  recruitment  dynamics  studies  alluded  to  the  concept  of  metapop¬ 
ulations,  population  structure  was  not  explicitly  considered  until  Roughgarden  and 
Iwasa  (1986)  and  Iwasa  and  Roughgarden  (1986).  Since  then,  there  has  been  much 
discussion  of  population  connectivity  and  the  importance  of  migration  for  structuring 
marine  populations. 

Despite  interest  in  population  connectivity  and  marine  metapopulations,  empiri¬ 
cal  studies  of  migration  rates  are  few  in  number  for  marine  invertebrates.  Much  of  the 
published  literature  focuses  on  migration  estimates  based  on  genetics  (e.g.  Grosberg 
and  Cunningham,  2001;  Gilg  and  Hilbish,  2003;  Sotka  and  Palumbi,  2006)  or  theoret- 
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ical  models  (e.g.  Incze  and  Naimie,  2000;  Roegner,  2000).  The  paucity  of  empirical 
migration  estimates  for  benthic  invertebrates  is  due  to  characteristics  of  the  stage 
in  which  dispersal  takes  place;  the  larval  stage  is  very  small  and  is  subject  to  high 
mortality  and  dilution  rates  (Thorson,  1950,  1966).  Despite  these  obstacles,  in  recent 
years  the  field  of  marine  ecology  has  sought  to  quantify  the  exchange  of  individuals 
among  geographically  separated  patches  in  a  metapopulation,  termed  larval  connec¬ 
tivity  (Cowen  et.  al..  2007).  In  this  thesis,  I  worked  towards  developing  the  use  of 
shell  chemistry  as  a  means  of  measuring  migration  rates,  and  therefore  connectivity, 
between  patches  of  the  softshell  clam,  Mya  arenana.  1  also  used  genetics  as  a  natural 
tag  of  environment  to  estimate  gene  flow.  Finally,  1  developed  and  implemented  a 
theoretical  model  to  explore  metapopulation  dynamics  and  understand  the  effects  of 
migration  on  metapopulation  growth  rate. 

1.2  Thesis  Goals 

Mv  thesis  began  as  method  development;  I  set  out  to  study  connectivity  and  retention 
of  the  softshell  clam  Mya  arenana  in  the  Nort  hwest  Atlantic  using  elemental  signa¬ 
tures  incorporated  into  the  larval  shell  as  tags  of  natal  habitat.  My  first  goal  was 
to  verify  that  the  instrument  commonly  used  in  other  natural  tagging  studies,  laser 
ablation  inductively  coupled  plasma  mass  spectrometry  (ICP-MS),  could  be  used  for 
analyzing  bivalve  larval  shell.  Results  from  my  experiment  (presented  in  Chapter 
2)  instead  indicated  that  this  instrument,  in  the  configuration  used,  was  not  able 
to  analyze  the  extremely  thin  first  larval  shell  of  M.  arenaria.  These  results  were 
not  expected,  and  I  was  forced  to  re-think  my  approach  to  my  thesis  questions.  My 
results  serve  to  discourage  ecologists  from  adopting  methods  that  have  proven  useful 
for  other  systems  without  first  carefully  considering  and  validating  the  method  for 
their  particular  species.  The  general  approach  of  using  biogenic  carbonate  chemistry 
for  connectivity  studies  is  still  useful;  different,  sensors  may  be  more  sensitive,  or  even 
a  different  application  of  ICP-MS  might  be  developed  to  analyze  the  surface  of  larval 
shells.  Although  I  chose  not  to  pursue  further  technique  development,  the  potential 
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for  bivalve  shell  chemistry  to  provide  ecologically  relevant  information — whether  for 
connectivity  or  other  environment  indicators — remains. 

If  biogenic  carbonate  is  to  be  used  successfully  as  an  environmental  indicator 
(for  connectivity  studies  or  otherwise),  it  is  essential  to  understand  the  effects  of 
environmental  and  biological  factors  on  elemental  composition.  In  Chapters  3  and 
4  report  on  experiments  I  conducted  to  understand  some  of  the  factors  potentially 
affecting  the  composition  of  M.  arenaria  shell.  Chapter  3  focuses  on  the  effects  of 
growth  rate  and  age  on  elemental  incorporation  into  shell,  and  Chapter  4  focuses  on 
temperature  and  salinity  effects  on  both  larval  and  juvenile  shell. 

When  I  discovered  that  I  would  not  be  able  to  answer  connectivity  questions 
about  M.  arenaria  using  laser  ablation  ICP-MS,  at  least  in  its  present  configuration, 
I  began  pursuing  population  genetics  as  a  different  type  of  natural  tag.  The  spatial 
and  temporal  scale  at  which  genetics  is  useful  for  parsing  out  population  variability 
depends  on  the  genetic  marker  chosen.  I  sequenced  the  mitochondrial  cytochrome 
oxidase  I  gene,  which  has  proven  useful  in  connectivity  studies  of  other  invertebrates. 
I  included  populations  in  the  northwest  Atlantic,  the  North  Sea,  and  the  Pacific  in  my 
study;  these  sites  represent  the  species’  native  and  introduced  range.  I  was  therefore 
able  to  assess  connectivity  on  larger  temporal  and  spatial  scales,  and  examine  patterns 
relating  to  M.  arenaria' s  natural  history.  These  results  are  reported  in  Chapter  5. 

The  variability  in  COI  that  I  was  able  to  detect  proved  informative  about  the 
history  of  M.  arenaria  over  global  spatial  scales.  Although  my  results  showed  that 
populations  in  New  England  were  not  isolated  from  one  another,  I  was  not  able  to 
address  larval  exchange  rates  over  ecologically  relevant  time  scales.  In  an  attempt  to 
return  to  my  original  questions  about  population  connectivity  and  metapopulation 
dynamics,  I  began  using  mathematical  models  of  M.  arenaria  metapopulations.  The 
matrix  population  model  presented  in  Chapter  6  began  by  my  adding  spatial  structure 
to  previous  work  on  M.  arenaria  by  Ripley  and  Caswell  (2006).  As  I  developed  the 
model,  it  metamorphosed  into  something  more  generally  applicable  to  any  metapop¬ 
ulation  subjected  to  stochastic  disturbance.  The  model  can  be  used  to  understand 
how  metapopulation  dynamics  are  affected  by  environmental  variability.  In  addition, 
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I  used  sensitivity  analysis  to  determine  the  effects  of  migration  on  metapopulation 
growth  rate.  I  found  that  since  migration  dictates  the  relative  distributions  of  indi¬ 
viduals  in  patches,  metapopulation  growth  rate  is  highly  sensitive  to  the  proportion 
of  individuals  moving  between  patches. 

Although  the  chapters  of  this  thesis  are  presented  in  chronological  order,  it  might 
be  more  useful  to  consider  the  results  of  Chapter  6  as  evidence  of  the  importance  of 
migration  for  determining  metapopulation  dynamics,  providing  justification  for  my 
attempts  to  study  dispersal  in  previous  data  chapters.  The  themes  and  questions  of 
this  thesis  are  generally  applicable  to  other  marine  benthic  invertebrates  and  lend 
knowledge  to  the  broader  concept  of  populat  ion  connectivity  and  the  techniques  used 
to  explore  this  important  ecological  concept  . 
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Laser  ablation  ICP-MS  analysis  of  larval  shell  in  softshell  clams 
(Mya  arenaria)  poses  challenges  for  natural  tag  studies 

C.A.  Strasser,  S.R.  Thorrold,  V.R.  Starczak,  and  L.S.  Mullineaux 

Woods  Hole  Oceanographic  Institution,  Woods  Hole,  MA  02543 

Abstract 

We  investigated  whether  laser  ablation  inductively  coupled  plasma  mass  spectrometry  (ICP-MS)  could  be  used 
to  quantify  larval  shell  compositions  of  softshell  clams,  Mya  arenaria.  The  composition  of  aragonitic  otoliths  has 
been  used  as  a  natural  tag  to  identify  natal  habitat  in  connectivity  studies  of  fish.  If  it  is  possible  to  measure  lar¬ 
val  shell  reliably,  this  technique  could  also  be  applied  to  marine  bivalves.  To  determine  whether  the  first  larval 
shell  (prodissoconch  I)  could  be  measured  independent  of  underlying  material,  we  conducted  laboratory  exper¬ 
iments  in  which  larval  M.  arenaria  were  exposed  to  enrichments  of  the  stable  isotope  UKBa  during  different  stages 
in  shell  development.  We  were  unable  to  isolate  the  chemical  signature  of  the  prodissoconch  I  from  subsequent 
life  stages  in  all  combinations  of  shell  preparation  and  instrument  settings.  Typical  instrument  settings  burned 
through  the  prodissoconch  I  on  a  post-settlement  juvenile  and  at  least  9  d  of  second  larval  shell  (prodissoconch  II) 
growth.  Our  results  suggest  instrumental  and  technical  improvements  are  needed  before  laser  ablation  ICP-MS 
can  be  useful  for  connectivity  studies  that  require  analysis  of  larval  shell  on  a  post-settlement  M.  arenaria  juve¬ 
nile.  Laser  burn-through  is  potentially  a  problem  in  any  connectivity  study  where  it  is  necessary  to  assay  the 
small  amounts  of  shell  material  that  are  deposited  before  a  larva  disperses  away  from  its  natal  location. 


Most  marine  benthic  invertebrate  life  cycles  include  a  plank¬ 
tonic  larval  phase  that  facilitates  dispersal  among  adult  popula¬ 
tions  (Thorson  1950).  Connectivity,  or  the  degree  to  which  geo¬ 
graphically  separated  populations  exchange  individuals,  is  an 
important  factor  in  the  spatial  population  dynamics  of  many 
marine  organisms  (Moilanen  and  Nieminen  2002).  An  under¬ 
standing  of  connectivity  in  marine  benthic  populations  is 
important  because  of  the  role  spatial  dynamics  play  in  fisheries 
management  and  the  design  and  implementation  of  marine 
protected  areas.  However,  studying  larval  dispersal  is  challeng¬ 
ing  due  to  small  larval  sizes,  high  dilution  rates,  and  high  larval 
mortality  rates  (Thorson  1950,  1966). 

In  recent  years,  the  use  of  artificial  and  natural  tags  to  track 
marine  larvae  has  been  explored  (e.g.,  Levin  1990;  Thorrold  et 
al.  2002).  One  type  of  natural  tag  that  may  be  useful  for  identi¬ 
fying  natal  origins  is  elemental  signatures  recorded  in  biogenic 
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carbonate.  This  technique  relies  on  the  observation  that  some 
elements  are  incorporated  into  the  calcium  carbonate  matrix 
in  amounts  that  are  related  to  the  dissolved  concentrations 
and  physical  properties  of  the  ambient  water  (e.g.,  Bath  et  al. 
2000;  Elsdon  and  Gillanders  2003;  Fowler  et  al.  1995;  Thorrold 
et  al.  1997;  Vander  Putten  et  al.  2000).  Provided  water  chem¬ 
istry  or  temperature  is  significantly  different  among  natal 
habitats,  such  variation  can  serve  as  a  natural  tag,  or  signature, 
of  the  geographic  origin  of  organisms.  The  use  of  geochemical 
signatures  in  fish  otoliths  as  natural  tags  for  population  stud¬ 
ies  is  well  established  (Campana  and  Thorrold  2001).  Recent 
efforts  have  expanded  the  use  of  elemental  tags  to  inverte¬ 
brates  including  decapods  (DiBacco  and  Levin  2000),  gas¬ 
tropods  (Zacherl  et  al.  2003),  bivalves  (Becker  et  al.  2005),  and 
cephalopods  (Arkhipkin  et  al.  2004). 

Most  studies  attempting  to  obtain  time-resolved  elemental 
signatures  from  calcified  tissues  have  used  laser  ablation  induc¬ 
tively  coupled  plasma  mass  spectrometry  (ICP-MS).  Conven¬ 
tional  solution-based  ICP-MS  analyses  are  generally  more  pre¬ 
cise  than  laser  ablation  assays  but  lack  the  ability  to  resolve 
elemental  signatures  from  individual  life  stages  (Campana  et  al. 
1997).  Laser  ablation  ICP-MS  is  particularly  useful  for  larval 
studies  because  it  allows  the  core  of  an  otolith  or  the  larval  shell 
of  a  juvenile  bivalve  to  be  targeted.  Yet,  surprisingly  few  studies 
have  empirically  tested  the  effective  spatial  resolution  of  laser 
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Prodissoconch  I  Prodissoconch  II  Juvenile  Shell 


Fig.  1.  Shell  development  in  M.  arenaria.  PI  and  Pll  larval  shells  are 
underlain  by  juvenile  and,  ultimately,  adult  shell. 

ablation  ICP-MS.  The  goal  of  this  study  was  to  determine  if 
techniques  that  have  proven  successful  in  otolith  studies  were 
also  useful  for  quantifying  larval  dispersal  of  the  softshell  clam, 
Mya  arenaria.  Specifically,  we  tested  if  it  was  possible  to  isolate 
and  measure  the  early  larval  shell  (prodissoconch  I)  of  M.  are¬ 
naria  using  laser  ablation  ICP-MS  without  contamination  from 
underlying  late  larval  and  juvenile  shell  material.  This  study  is 
an  important  step  in  developing  the  use  of  elemental  signatures 
in  mollusc  shells  as  tags  of  natal  habitat  because  isolation  of  the 
early  larval  shell  is  critical  for  connectivity  studies. 

Materials  and  procedures 

Study  species — Adult  M.  arenaria  broadcast  spawn  their 
gametes  into  the  water,  where  fertilized  eggs  develop  into 
non-feeding  shell-less  trochophore  larvae.  Within  24  h  of 
spawning,  the  trochophore  metamorphoses  into  a  veliger  and 
lays  down  its  first  larval  shell,  the  prodissoconch  1  (hereafter  PI), 
which  is  approximately  90  urn  in  diameter  (Loosanoff  et  al. 
1966).  After  forming  the  PI,  the  veliger  begins  feeding  in  the 
plankton.  Over  the  next  1  to  3  weeks,  it  lays  down  the  second 
larval  shell,  the  prodissoconch  II  (hereafter  PII),  which  grows 
to  diameters  of  150  to  200  pm  (Loosanoff  et  al.  1966).  The 
veliger  then  finds  suitable  habitat,  settles  to  the  bottom,  and 
metamorphoses  into  a  juvenile  M.  arenaria.  Juvenile  and, 
eventually,  adult  shell  is  laid  underneath  larval  shell,  extend¬ 
ing  away  from  the  umbo  (Fig.  1).  Larval  shells  are  retained  dur¬ 
ing  clam  growth  and  development,  and  the  PI  can  be  identi¬ 
fied  easily  on  newly  settled  juveniles.  Estimated  water 
retention  times  for  typical  New  England  estuaries  are  2  days  to 
a  week  (Asselin  and  Spaulding  1993;  Brooks  et  al.  1999;  Shel¬ 
don  and  Alber  2002).  If  larvae  are  essentially  passive  particles, 
then  they  are  also  retained  in  estuaries  for  2  days  to  a  week 
after  spawning,  which  would  be  sufficient  for  the  PI  larval 
shell  to  form  before  potential  dispersal  occurs.  If  the  composi¬ 
tion  of  the  larval  shell  on  a  juvenile  M.  arenaria  can  be  mea¬ 
sured  accurately,  we  may  be  able  to  answer  questions  about 
larval  retention  and  dispersal. 

Claw  rearing — M.  arenaria  veligers  were  obtained  from  the 
Eastham  Hatchery  Facility  in  Eastham,  Massachusetts,  within 
6  h  of  spawning  and  transported  to  Woods  Hole,  Massachusetts. 


Unspiked 

No  shell  spiked 

Early  Spike 

Pll,  juvenile  shell 
spiked 

Late  Spike 

Some  juvenile  shell 
spiked 

Fig.  2.  Location  of  the  ,i8Ba  spike  treatments  relative  to  shell  develop¬ 
ment.  Shaded  regions  show  areas  of  spiked  shell. 

Approximately  400, (XX)  individuals  were  divided  into  nine 
identical  covered  10  L  tanks.  Tanks  were  filled  with  20°C  fil¬ 
tered  seawater  (salinity  -  32).  Assuming  roughly  10%  mortal¬ 
ity  during  transport  from  Eastham  to  Woods  Hole,  initial  lar¬ 
val  counts  were  approximately  40, (XX)  per  tank.  Each  tank  was 
randomly  assigned  one  of  three  treatments  described  below, 
with  three  replicate  tanks  per  treatment. 

Several  larvae  were  subsampled  from  each  tank  after  48  h, 
and  the  PI  was  present  in  all  sampled  larvae.  We  then  exposed 
clams  to  one  of  three  treatment  conditions  that  differed  in  the 
timing  of  exposure  to  the  stable  isotope  ,  ,MBa  (Fig.  2).  Tanks  in 
the  first  treatment  received  no  I,KBa  spike  and  therefore  the 
water  had  a  natural  MKBa:,<7Ba  ratio  ("Unspiked").  Tanks  in  the 
second  treatment  received  a  HMBa  spike  of  25  pg  fr1  at  day  2, 
immediately  after  the  PI  had  set  ("Early  Spike").  Tanks  in  the 
third  treatment  received  a  MKBa  spike  of  25  pg  g-1  at  day  9,  7  d 
after  the  PI  set  ("Late  Spike").  By  day  9,  most  of  the  larvae 
appeared  to  have  set  the  Pll  and  become  competent  to  settle 
based  on  their  shell  size,  tendency  to  accumulate  at  the  bottom 
of  the  tank,  and  use  of  their  foot  for  locomotion.  This  is  faster 
than  the  estimated  two  to  three  weeks  residence  time  of  M.  are¬ 
naria  veligers  in  the  field.  Laboratory  clams  typically  grow  and 
mature  faster  compared  with  those  in  the  field  because  of  opti¬ 
mized  rearing  conditions  (H.  Lind,  pers.  comm.). 

Larvae  swam  freely  in  tanks  for  the  first  10  d,  after  which 
time  they  were  contained  in  mesh  sieves  suspended  in  the 
tanks.  Water  pumps  were  used  to  create  a  downward  current 
that  encouraged  larvae  to  settle  out  of  the  water  column.  Sub¬ 
sampling  of  tanks  indicated  that  larvae  metamorphosed 
within  2  to  3  d  after  downwelling  began.  Clams  were  fed 
Isochrvsis  sp.  and  Pavlova  sp.  algae  at  least  once  every  2  d,  and 
the  water  in  each  tank  was  changed  every  2  d.  In  the  Early  and 
Late  Spike  tanks,  a  MMBa  spike  concentration  of  25  pg.fr1  was 
maintained  for  the  duration  of  the  experiment  by  adding  1  ,HBa 
immediately  after  each  water  change.  After  60  d,  clams  ranged 
in  size  from  0.5  mm  to  5  mm.  They  were  removed  from  their 
tanks  and  frozen  until  prepared  for  ICP-MS  analysis.  Mortality 
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rates  were  high  in  all  tanks,  and  seven  of  the  nine  tanks  had 
survivors:  three  replicate  tanks  from  the  Unspiked  treatment 
(//  =35 ,  n  -  33,  n  -  7),  two  replicate  tanks  from  Early  Spike 
treatment  (n  =  75,  n  =  22),  and  two  replicate  tanks  from  the 
Late  Spike  treatment  ( n  =  14,  n  =  8). 

Shell  preparation — Shells  were  cleaned  thoroughly  using 
techniques  developed  for  foraminiferan  tests  (Boyle  1981) 
with  modifications  specifically  for  M.  arenaria.  Clams  were 
placed  in  individual  acid-washed  vials  using  acid-washed  plas¬ 
tic  forceps,  and  sonicated  briefly  to  remove  tissue  and  debris. 
Individuals  were  rinsed  three  times  with  ultrapure  H20,  then 
soaked  for  10  min  at  80°C  in  1%  H202  solution  buffered  with 
suprapur  NaOH  (1  N)  to  remove  organic  material.  Afterward, 
shells  were  rinsed  three  times  with  ultrapure  H,0,  transferred 
to  clean,  acid-washed  vials,  rinsed  four  times  with  ultrapure 
H20,  then  left  to  dry  overnight  under  a  laminar  flow  hood. 
Shells  were  inspected  after  cleaning  and  drying  for  any 
remaining  organic  material;  dried  tissue  was  removed  in  the 
few  cases  where  it  was  present. 

Shells  were  mounted  for  laser  ablation  ICP-MS  analysis  on 
glass  slides  using  Devcon©  Super  Glue.  One  valve  of  each  clam 
was  oriented  so  that  the  retained  larval  shell  was  accessible  by 
the  laser  (i.e.,  umbo  facing  upward),  while  the  second  valve 
was  oriented  concave  side  down  so  that  the  laser  could  access 
the  flattest  portion  of  the  juvenile  shell.  All  cleaning  and 
preparations  were  conducted  in  a  Class  100  clean  room. 

Laser  ablation  ICP-MS  analysis — Shell  material  was  analyzed 
using  a  Thermo-Electron  Element2  ICP-MS  coupled  to  a  New 
Wave  Research  213  nm  laser  ablation  system.  Vaporized  mate¬ 
rial  from  the  ablation  was  transported  via  a  helium  gas  stream 
to  the  dual-inlet  quartz  spray  chamber  where  it  was  mixed 
with  a  2%  HNO,  aerosol  from  a  self-aspirating  PFA  20  pL  min  ' 
nebulizer.  The  analyte  was  then  transported  to  the  ICP-MS  via 
an  argon  carrier  gas.  An  aragonite  otolith  reference  material 
dissolved  in  2%  HNO,  and  diluted  to  a  Ca  concentration  of 
40  ng  g-'  (Yoshinaga  et  al.  2000)  and  a  2%  HNO^  blank  solu¬ 
tion  were  run  periodically  to  correct  for  mass  bias  drift.  We 
measured  4KCa,  l  lKBa,  and  ,;<7Ba  in  low  resolution  mode.  Iso¬ 
tope  ratios  of  UKBa  to  1 17Ba  were  calculated  using  mass  bias  cor¬ 
rection  calculated  from  calibration  standards  that  we  assumed 
contained  natural  °8Ba:u7Ba  ratios.  Blank  intensity  averages 
and  standard  deviations  were  calculated  for  each  of  the  three 
analysis  periods  (Oct  2005,  Jan  2006,  Feb  2006).  Limits  of 
detection  (LOD)  were  calculated  as  the  ratio  of  three  standard 
deviations  of  the  blank  intensity  to  the  average  blank-sub¬ 
tracted  sample  intensity.  For  the  three  time  periods  over 
which  samples  were  analyzed,  sample  intensities  were  >  100, 

>  690,  and  >  2800  times  the  detection  limit  for  4KCa;  >  100, 

>  30,  and  >  20  times  the  detection  limit  for  ,17Ba;  and  >  290, 

>  127,  and  >  30  times  the  detection  limit  for  1  <KBa,  respectively. 

Measurements  were  taken  from  two  positions  on  each 
clam:  one  on  the  larval  shell  at  the  umbo,  and  a  second  on 
juvenile  shell  material  600  fim  away  from  the  umbo  (umbo 
and  juvenile  measurements,  respectively).  Shell  composition 


Table  1.  Laser  ablation  ICP-MS  settings  used  to  analyze 
M.  arenaria  shell 


Raster  setting 

Spot  setting 

Laser  pattern 

70  nm2  raster 

80  nm  spot 

Laser  power 

80% 

50% 

Beam  size 

25  nm 

80  nm 

Repetition  rate 

10  Hz 

5  Hz 

Scan  speed 

12  nm  r' 

Line  spacing 

10  nm 

for  each  clam  was  analyzed  using  one  of  two  laser  settings 
(Table  1).  In  the  first  case  (Raster  Setting),  the  laser  removed 
shell  material  systematically  from  a  70  pm  x  70  pm  square 
raster  pattern.  Instrument  and  laser  settings  were  chosen  to 
maintain  4HCa  counts  above  two  million  counts  s-'.  For  the 
second  setting  (Spot  Setting),  the  laser  was  focused  on  a  single 
spot  with  a  beam  diameter  of  80  pm,  and  laser  energy  was 
reduced  to  50%  power.  The  laser  pulsed  on  the  spot  repeat¬ 
edly,  boring  down  into  shell  material.  Therefore,  at  least  the 
first  few  seconds  of  measurement  at  the  umbo  should  reflect 
larval  shell;  afterward  the  laser  was  expected  to  penetrate 
underlying  material.  This  setting  was  used  to  explore  whether 
the  composition  of  thin  larval  shell  could  be  quantified  at 
reduced  laser  power  with  minimal  material  removed  (i.e.,  only 
a  few  seconds  of  material  ablation). 

Rationale  for  experimental  setup — Umbo  and  juvenile  shell 
were  each  measured  on  clams  from  the  three  Unspiked  treat¬ 
ment  tanks  (/;  =  13,  n  =  9,  n  =  6)  to  test  whether  uptake  dif¬ 
fered  between  umbo  and  juvenile  areas  of  the  shell.  No  statis¬ 
tically  significant  differences  were  found  for  the  three  tanks 
combined  (paired  t  test,  t  =  0.83,  //  =  28 ,  P  -  0.42)  so  only 
umbo  measurements  were  used  to  simplify  analysis. 

Shell  ratios  of  1,MBa:°7Ba  were  expected  to  fall  into  one  of 
three  categories:  a  natural  ratio,  a  fully  UKBa-spiked  elevated 
ratio,  or  somewhere  between  these  two  end  members.  The 
natural  ratio  was  determined  by  measuring  shell  that  was  not 
spiked  (i.e.,  clams  from  the  Unspiked  treatment).  The  fully 
spiked  ratio  was  determined  for  each  clam  individually  in  the 
Early  and  Late  Spike  treatments  by  measuring  juvenile  shell. 
Umbo  measurements  of  clams  from  Early  Spike  and  Late  Spike 
treatments  were  expected  to  show  one  of  three  results, 
depending  on  the  depth  of  laser  penetration  (Fig.  3):  (A)  no 
laser  burn-through:  umbo  measurements  from  Early  and  Late 
Spike  treatments  did  not  significantly  differ  from  unspiked 
shell;  (B)  some  laser  burn-through:  umbo  measurements  from 
the  Early  Spike  treatment  have  a  significantly  higher  ratio 
than  unspiked  shell;  umbo  measurements  from  the  Late  Spike 
treatment  have  a  ratio  that  does  not  differ  significantly  from 
unspiked  shell;  (C)  extensive  laser  burn-through:  umbo  mea¬ 
surements  from  Early  and  Late  Spike  treatment  have  ratios 
that  do  not  differ  significantly  from  spiked  juvenile  shell. 

If  laser  burn-through  does  occur,  we  may  be  able  to  cor¬ 
rect  for  it  by  subtracting  out  the  signal  from  any  underlying 
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the  umbo  measurement.  Equation  1  was  rearranged  to  deter¬ 
mine  the  proportion  of  spiked  material  for  each  clam: 


A  No  laser  bum-through 


v 


£  20 


Early  Spike  Late  Spike 
umbo  umbo 


B  Some  laser  bum-through 


Early  Spike  Late  Spike 
umbo  umbo 

C  Extensive  laser  bum-through 


J 


Early  Spike  Late  Spike 
umbo  umbo 


Fig.  3.  Three  possible  scenarios  for  extent  of  laser  burn-through.  Graphs 
are  the  expected  '  '"Ba: 1 ,7Ba  ratios  for  umbo  measurements  from  the  three 
treatments,  depending  on  depth  of  laser  penetration:  (A)  no  laser  bum- 
through,  (B)  some  laser  burn-through,  and  (C)  extensive  laser  burn- 
through.  Shaded  regions  of  the  shell  graphics  indicate  spiked  material. 
The  shaded  horizontal  line  in  each  graph  at  Expected  mBa:,J7Ba  -20 
shows  fully  spiked  shell,  while  the  line  at  Expected  mBa:n;Ba  -6  shows 
the  natural  n8Ba:H7Ba  ratio. 


shell  that  might  be  contaminating  our  results.  To  determine 
if  this  was  possible,  we  used  a  simple  mixing  equation  (Eq.  1) 
to  calculate  the  proportion  of  spiked  and  unspiked  material 
ablated  for  each  clam's  umbo  measurement  from  the  spiked 
treatments: 


(1) 


where  X  is  the  proportion  of  spiked  versus  unspiked  material 
ablated  from  the  umbo  of  the  clam,  (,:<MBa:M7Ba)Spjked  is  the 
ratio  from  the  juvenile  measurement  of  the  clam, 
(U8Ba:M7Ba)Unsplkw)  is  the  average  ratio  from  umbos  of  clams  in 
the  Unspiked  treatment,  and  (MKBa:u7Ba)Umbii  is  the  ratio  from 
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Editing  profiles  from  Spot  Setting — The  Spot  Setting  was  used 
to  determine  if  unspiked  material  could  be  detected  by  pro¬ 
ducing  a  Ba  isotope  depth  profile  centered  on  the  PI  position 
of  the  shell.  Data  were  initially  acquired  from  a  2%  HNOt 
blank  solution,  with  laser  ablation  beginning  approximately 
18  s  into  data  collection.  The  resulting  intensity  profiles  were 
edited  as  follows  to  determine  blank  intensities  and  shell 
material  intensities  for  the  three  isotopes  of  interest.  Blank 
values  were  calculated  as  the  average  of  the  first  18  s  of  analy¬ 
sis.  We  started  collecting  data  for  shell  material  once  a  data 
point  was  SOM)  higher  in  4KCa  intensity  than  its  predecessor  for 
a  total  of  20  s.  Detection  limits  were  then  calculated  for  the 
three  isotopes  based  on  the  blank  intensities,  and  only  data 
with  intensities  at  least  20%  above  the  detection  limit  for  all 
three  isotopes  were  included  in  analyses  (Fig.  4).  The  criterion 
used  to  remove  data  below  the  detection  limit  resulted  in 
some  profile  sequences  being  shorter  than  others;  time  pro¬ 
files  ranged  from  9  to  20  data  points. 


Assessment 

Differences  between  treatments — Using  the  Raster  Setting,  the 
average  (±  SD)  ,3MBa:,37Ba  ratio  in  shell  material  was  6.4  ±  1.0  for 
Unspiked  umbo  (//=  19),  20.3  ±  6.1  for  Early  Spike  umbo  (//=  12), 
and  15.5  ±  2.3  for  Late  Spike  umbo  (//  =  11)  (Fig.  5).  Using  the 
Spot  Setting,  the  average  (±  SD)  ratio  was  5.8  ±  1.1  for  unspiked 
umbo  (//=  8),  20.0  ±  1.9  for  Early  Spike  umbo  (//=  14),  and  15.5 
±  3.9  for  Late  Spike  umbo  (//  =  3).  To  determine  if  laser  burn- 
through  occurred,  we  first  tested  whether  there  was  significant 
variation  in  mean  umbo  measurements  among  treatments  and 
among  tanks  within  treatments  using  a  mixed  model  nested 
ANOVA  (Table  2).  Under  both  Raster  and  Spot  Settings,  mean 
umbo  measurements  were  significantly  different  between  the 
three  treatments  (Raster,  fzx,  =  56.4,  P<  0.(XX)1 ;  Spot,  ,,  =  72.5, 
P  <  0.0001),  with  Unspiked  ratio  <  Late  Spike  <  Karly  Spike 
(Tukey-Kramer  Test).  We  found  no  significant  differences  among 
tanks  within  treatments  for  the  raster  setting  (/>=  0.88;  number 
of  clams  per  replicate  tank:  n  =7,  //  =  6,  //  =  6  for  Unspiked;  n  -  6, 
//  =  6  for  Early  Spike;  //  =  6,  //  =  5  for  Late  Spike).  We  were  not 
able  to  test  for  differences  among  tanks  within  treatments  for  the 
spot  setting  due  to  insufficient  numbers  (number  of  samples;  n 
=  5,  //  =  3  for  Unspiked;  //  =  14  for  Early  Spike;  //  =  3  for  Late 
Spike).  Hereafter  measurements  from  tanks  of  the  same  treat¬ 
ment  were  pooled.  Because  both  Early  and  Late  Spike  umbo 
ratios  differed  significantly  from  Unspiked  umbo,  we  rejected 
laser  burn-through  scenarios  A  and  B  and  concluded  that  exten¬ 
sive  laser  burn-through  was  occurring  (Fig.  3C). 
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Fig.  4.  Representative  unedited  intensity  profiles  for  individual  clams  from  laser  ablation-ICP-MS  measurements  using  the  Spot  Setting  for  (A)  48Ca  and  (B) 
H8Ba.  The  data  in  the  first  shaded  region  were  used  to  calculate  the  blank  levels,  and  the  data  in  the  second  shaded  region  were  excluded  from  analysis. 


Ratios  for  individual  clams — The  extent  of  laser  burn-through 
was  addressed  by  comparing  umbo  measurements  with  spiked 
juvenile  shell  for  each  clam  using  a  paired  t  test.  If  the  umbo 
ratio  was  not  significantly  different  from  spiked  juvenile  shell, 
it  suggests  that  a  measurable  quantity  of  unspiked  shell  was  not 
present  at  the  umbo.  For  the  Raster  Setting,  umbo  mea¬ 
surements  were  not  significantly  different  from  juvenile  mea¬ 
surements  in  Unspiked  clams  (/=  -0.45,  n  =  19,  P-  0.66)  or 
Early  Spike  clams  (/  =  -1.0,  n  =  12,  P  -  0.34).  In  Late  Spike 
clams,  the  umbo  ,<8Ba:u7Ba  ratio  of  15.6  was  slightly  lower 
than  the  juvenile  shell  ratio  of  18.3  ( t-  2.7,  ti-  11,  P<  0.05). 
Results  were  similar  using  the  Spot  Setting.  Unspiked  and  Early 
Spike  umbo  means  were  not  significantly  different  from  juve¬ 
nile  measurements  (/=  1.6,  n  -  9,  P=  0.15;  and  t  -  0.85,  n  =  14, 
P-  0.41,  respectively),  whereas  Late  Spike  mean  umbo  and 
juvenile  ratios  differed,  although  sample  size  was  small  (/=-7.0, 
n  =  3,  P<  0.05).  The  umbo  and  juvenile  ratios  are  plotted  for 
each  clam  from  the  Raster  Setting  (Fig.  6).  Results  were  similar 
for  Spot  Setting  analyses  and  are  not  presented.  The  data  were 
again  consistent  with  a  scenario  of  extensive  laser  bum- 
through  (Fig.  3C).  Umbo  shell  ratios  tended  to  be  more  similar 
to  juvenile  shell  ratios  than  to  unspiked  shell,  indicating  that 
umbo  measurements  were  not  isolating  the  unspiked  larval 
shell  from  underlying  shell  material. 

For  some  clams,  umbo  ratios  were  unexpectedly  higher 
than  juvenile  ratios  (Fig.  6).  In  theory,  all  juvenile  clams  expe¬ 
rienced  the  maximum  amount  of  spike,  but  ratios  in  juvenile 
shell  were  variable  among  clams,  ranging  from  9.4  to  22.3. 
Some  of  this  variability  might  be  attributed  to  physiological 
differences  among  clams,  but  much  of  it  likely  originates  from 
variability  in  the  laser  ablation  ICP-MS  measurements.  For 
instances  when  umbo  ratios  were  higher  than  juvenile  ratios 
for  the  same  clam,  the  magnitude  of  the  difference  was  within 
the  range  of  variability  expected  based  on  the  range  seen  in 
juvenile  shell  ratios. 


Late  Spike  umbo  ratios  were  slightly  lower  than  Early  Spike 
umbo  ratios,  indicating  that  a  larger  proportion  of  unspiked 
shell  was  being  ablated  in  Late  Spike  clams.  Clams  in  this  treat¬ 
ment  were  allowed  to  set  shell  for  9  d  without  spike;  however 
measurements  from  the  umbo  still  showed  only  a  small  differ¬ 
ence  from  fully  spiked  shell.  This  is  a  surprising  result  with 
implications  for  the  usefulness  of  ICP-MS  in  connectivity  stud¬ 
ies.  Measurements  taken  at  the  umbo  may  not  represent  natal 
habitat  signatures  alone,  but  rather  the  integration  of  habitats 
encountered  by  the  organism  over  more  than  9  d  of  life. 

Time  profiles  of  umbo  measurements — Spot  Setting  data  were 
used  to  investigate  whether  unspiked  larval  shell  could  be 


Fig.  5.  Average  umbo  1i8Ba:li7Ba  ratios  (±  standard  error)  for  clams  mea¬ 
sured  using  Raster  and  Spot  Settings  in  the  three  treatments.  The  shaded 
horizontal  line  shows  the  natural  ,J8Ba:137Ba  ratio. 
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Table  2.  Results  from  ANOVA  testing  for  differences  among 
treatment  mean  ratios 


Source 

df 

MS 

F 

_ P 

Raster  Setting 

Treatment 

2 

760 

56.4 

<  0.0001 

Tank  (Treatment) 

3 

5.36 

0.397 

0.877 

Error 

35 

13.5 

Spot  Setting 

Treatment 

2 

497 

72.5 

<  0.0001 

Error 

23 

6.86 

For  the  raster  setting,  results  are  from  a  mixed  model  nested  ANOVA  test¬ 
ing  for  differences  among  treatment  means  and  for  variation  among  tanks 
nested  within  treatments.  For  the  spot  setting,  results  are  from  a  one-way 
ANOVA;  data  from  tanks  were  combined  to  have  sufficient  sample  size  to 
test  for  differences  among  treatment  means  (see  Methods). 

detected  in  the  first  seconds  of  laser  ablation.  This  would  be 
evident  if  initial  MKBa:M7Ba  ratios  were  near  natural  levels,  or 
at  least  significantly  lower  than  subsequent  ratios.  Visual 
inspection  of  depth  profiles  suggested  that  unspiked  shell  was 
not  being  detected  early  during  ablation  (Fig.  7).  To  statisti¬ 
cally  test  this,  we  averaged  the  first  five  points  of  shell  ablation 
data  and  compared  the  value  with  an  average  of  the  remain¬ 
ing  15  points  using  a  paired  /  test.  Five  points  were  averaged 
to  represent  initial  shell  because  single  data  points  are  highly 
variable  and  subject  to  small  fluctuations  in  the  amount  of 
material  ablated  and  instrument  sensitivity  (Guillong  et  al. 
2001;  Russo  1995).  If  a  significant  proportion  of  material 
ablated  early  was  composed  of  unspiked  shell,  the  first  five 
points  should  be  significantly  lower  than  the  last  15.  In  all 
treatments,  the  first  five  points  were  not  significantly  lower 
than  the  last  15  (Unspiked:  /=  1.2,  n  =  7,  P-  0.29;  Early  Spike: 
/=  1.979,  //=  13,  P=  0.071;  Late  Spike:  /=- 1.2,  //=  5,  />=.30). 
The  lack  of  unspiked  material  might  be  due  to  the  over¬ 
whelming  signal  of  underlying  spiked  shell,  or  that  the  points 
representing  unspiked  shell  fell  below  the  >  20%  detection 
limit  cutoff  and  were  excluded.  In  either  case,  our  results  sug¬ 
gested  that  unspiked  larval  shell  was  not  detectable  even  when 
laser  power  was  reduced  and  minimal  material  ablated. 

Neither  Early  nor  Late  Spike  treatments  were  spiked  during 
PI  formation.  As  a  result,  initial  shell  measurements  at  the 
umbo  should  have  reflected  unspiked  shell  if  only  the  PI  was 
measured.  Measurements  taken  later  (deeper)  were  expected  to 
have  an  increase  in  the  ratio  of  ,<KBa:u7Ba,  and  the  timing  of 
that  increase  would  be  related  to  the  time  when  the  spike  was 
added  during  shell  formation.  The  absence  of  detectable 
unspiked  shell  at  the  beginning  of  ablation  for  both  Early  and 
Late  Spike  treatments  suggests  that  laser  ablation  does  not  iso¬ 
late  and  measure  the  first  larval  shell  accurately.  Furthermore, 
based  on  results  from  Late  Spike  clams,  the  first  few  seconds  of 
laser  ablation  removes  at  least  9  d  worth  of  shell  growth  from 
laboratory-reared  clams. 

Although  we  were  unsuccessful  at  measuring  the  larval  shell 
in  our  study,  there  is  no  evidence  that  the  PI  larval  shell  is  meta- 


Flg.  6.  U8Ba:n/Ba  ratios  for  each  clam,  for  each  of  the  three  treatments 
as  measured  by  the  Raster  setting:  (A)  Unspiked,  (B)  Early  Spike,  and  (C) 
Late  Spike.  The  shaded  horizontal  line  in  each  graph  shows  the  natural 
U8Ba:li7Ba  ratio 


bolically  reworked,  absent  from  the  juvenile  shell,  or  otherwise 
unusable  for  connectivity  studies.  However,  our  study  demon¬ 
strated  that  regardless  of  the  presence,  absence,  or  inert  qualities 
of  the  larval  shell,  we  were  unable  to  detect  or  measure  the  lar¬ 
val  shell  with  the  settings  and  instrumentation  reported. 

Proportion  spiked  material — We  calculated  the  proportion  of 
spiked  material  for  Early  and  Late  Spike  umbo  measurements 
from  the  Raster  Setting  using  Eq.  2  (Fig.  8).  Average  propor¬ 
tions  of  spiked  material  (±  SD)  were  1.3  ±  0.7  (//=  12)  for  Early 
Spike  umbo  measurements  and  0.9  ±  0.5  (//=  11)  for  Late  Spike 
umbo.  A  one-sample  t  test  showed  no  significant  difference 
between  mean  proportions  and  unspiked  material,  i.e.,  the 
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Fig.  7.  Three  representative  edited  time  profiles  of  138Ba:137Ba  for  individual  clams  using  the  Spot  Setting:  (A)  Early  Spike  umbo,  (B)  Late  Spike  umbo, 
(C)  spiked  juvenile  shell,  and  (D)  Unspiked  umbo.  The  horizontal  gray  line  indicates  the  natural  n8Ba:n;Ba  ratio.  Profiles  not  shown  had  similar  patterns. 


mean  proportion  did  not  differ  from  1.0  (Early  Spike:  t-  1.5, 
P-  0.16;  Late  Spike:  t=  -0.7,  />=.48).  Results  were  similar  for 
Spot  Setting  clams  and  are  not  shown. 

Variability  in  the  ratios  of  umbo  measurements  among  indi¬ 
viduals  from  the  two  spiked  treatments  resulted  in  a  high  degree 
of  variability  in  the  calculated  proportion  of  spiked  material.  We 
were  able  to  detect  this  variability  and  quantify  the  amount  of 
spiked  versus  unspiked  shell  since  material  underneath  the  lar¬ 
val  shell  was  tagged  using  UHBa.  However,  this  variability  is  not 
easily  predicted  or  quantified  in  field  samples  of  M.  arenaria 
since  there  is  no  known  element  with  different,  and  constant, 
concentrations  in  larval  versus  juvenile  shell.  Therefore,  we 
concluded  that  correcting  for  laser  burn-through  using  Eq.  2 
likely  would  not  be  useful  for  field  samples. 

Discussion 

Implications  for  connectivity— 'We  were  unable  to  isolate  and 
measure  the  larval  shell  of  M.  arenaria  using  laser  ablation 
1CP-MS  settings  designed  for  minimal  shell  removal.  The  laser 
consistently  burned  through  the  larval  shell  and  into  under¬ 
lying  late-stage  larval  and  juvenile  material.  In  addition,  the 
proportion  of  underlying  shell  ablated  was  too  variable  to 
allow  for  any  reliable  correction  to  account  for  the  burn- 
through.  Our  study  suggests  laser  burn-through  is  a  significant 
problem  for  connectivity  studies  of  M.  arenaria ,  and  possibly 
of  other  molluscan  species  that  may  spend  only  a  short  time 
in  their  natal  location. 


New  England  estuaries  tend  to  have  residence  times  on  the 
order  of  2  d  to  more  than  a  week  (Asselin  and  Spaulding  1993; 
Sheldon  and  Alber  2002),  suggesting  that  larvae  may  experi¬ 
ence  their  natal  habitat  for  as  few  as  48  h.  Our  results  indicate 
that  we  burned  through  at  least  9  d  of  shell  growth  in  the 
umbo  region  of  laboratory-reared  juvenile  M.  arenaria.  The 
number  of  days  of  growth  ablated  depends  on  shell  growth 
rate,  which  varies  depending  on  field  environmental  condi¬ 
tions  and  laboratory  rearing  conditions  (La Valley  2001). 
Although  there  may  be  some  differences  between  growth  rates 
in  the  field  and  those  of  clams  in  our  experiment,  our  results 
show  that  laser  ablation  removes  significantly  more  shell  than 
what  is  laid  in  the  first  24  to  48  h,  when  clams  are  most  likely 
retained  in  their  natal  habitat. 

Laser  burn-through  appears  to  be  a  substantial  problem  for 
connectivity  studies  of  organisms  with  rapid  dispersal,  but  it 
may  be  less  problematic  in  other  scenarios.  Species  that  brood 
their  larvae  may  produce  more  shell  material  before  potential 
dispersal,  making  isolation  of  natal  habitat  signatures  possible. 
Similarly,  larvae  that  are  broadcast  spawned  but  retained  in 
their  natal  habitat  for  most  of  their  larval  stage  due  to  physics 
might  have  time  to  deposit  sufficient  shell  for  reliable  mea¬ 
surement.  However,  requiring  that  organisms  experience  their 
natal  habitat  for  long  periods  of  time  limits  the  species  we  can 
study  using  natural  tagging  techniques. 

Alternative  approaches — In  this  study,  we  chose  to  explore  laser 
bum-through  issues  relating  to  laser  ablation  ICP-MS  because  it 
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Fig.  8.  Proportion  spiked  material  for  umbo  measurements  of 
M8Ba-treated  clams  using  the  Raster  Setting,  for  (A)  Early  Spike  and 
(B)  Late  Spike  treatments. 


is  an  instrument  commonly  used  for  carbonate  analysis  and  has 
been  chosen  by  other  molluscan  tagging  studies.  Our  results  sug¬ 
gest  that  other  techniques  and  instrument  configurations  should 
be  tested  to  determine  whether  they  might  more  reliably  isolate 
the  thin  larval  components  of  shell.  For  instance,  researchers 
have  reported  that  excimer  lasers  operating  at  193  nm  generate 
shallower  craters  than  we  were  able  to  achieve  in  the  present 
study  (e.g.,  Patterson  et  al.  2005),  while  the  burn-through  prob¬ 
lem  is  likely  more  pronounced  with  instalments  using  lasers 
with  longer  wavelengths  (e.g.,  FitzGerald  et  al.  2004;  Jones  and 
Chen  2003).  We  are  unaware  of  any  formal  studies  comparing 
crater  depths  in  carbonates  for  lasers  of  different  wavelengths. 
Gonzalez  et  al.  (2002)  showed,  however,  that  crater  shapes  and 
depths  produced  in  glass  standards  were  similar  for  193  nm  and 
213  nm  lasers,  and  both  of  these  lasers  produced  shallower  and 
more  regular  craters  than  266  nm  lasers. 


Instrument  sensitivity  also  plays  a  significant  role  in  deter¬ 
mining  the  amount  of  material  that  needs  to  be  transported  to 
the  plasma  to  make  a  sufficiently  precise  and  accurate  mea¬ 
surement  of  shell  chemistry.  We  did  not  experiment  with  run¬ 
ning  the  ICP-MS  in  dry  plasma  mode  with  a  desolvating  neb¬ 
ulizer,  because  in  our  experience  such  changes  lead  to  only 
small  increases  in  sensitivity.  Barnes  et  al.  (2004)  reported 
that  combining  a  laser  with  an  ICP  ionization  source  and  a 
Mattuch-Herzog  mass  spectrograph  can  dramatically  increase 
limits  of  detection  (i.e.,  instrument  sensitivity)  and  allow  for 
accurate  measurement  of  material  from  single  laser  pulses 
(5  nm  depth  per  pulse).  However,  development  of  such 
unique,  customized  systems  designed  for  specific  types  of 
analysis  is  difficult  and  expensive,  making  it  an  intractable 
solution  for  most  researchers. 

Other  instalments  that  may  warrant  further  investigation 
include  secondary  ionization  mass  spec  trometry  (SIMS)  and  pro¬ 
ton  induced  x-ray  emission  (PIXE).  Although  they  have  been 
used  successfully  in  otolith  studies,  LOD  and  penetration  depths 
of  the  secondary  ion  beam  (SIMS)  and  x-rays  (PIXE)  suggest  that 
neither  instalment  is  likely  to  isolate  larval  shell  material  effec¬ 
tively  (Campana  et  al.  1997).  Another  potential  problem  with 
SIMS  and  PIXE  is  that  present  configurations  are  not  able  to  ana¬ 
lyze  trace  elements  such  as  Ba,  Cd,  Mn,  and  Pb  as  effectively  as 
ICP-MS  (Campana  1999).  Such  elements  have  proven  useful  for 
distinguishing  different  habitats  in  natural  tagging  studies  of 
fishes  (e.g.,  Patterson  et  al.  2005;  Thorrold  et  al.  1998;  Thorrold 
et  al.  2001).  Finally,  both  SIMS  and  PIXE  require  the  surface  of 
the  analyte  to  be  flat.  For  bivalve  larval  shell  analysis,  this 
requires  that  the  juvenile  shell  be  embedded  in  epoxy,  then 
cross-sectioned  so  that  the  larval  shell  is  visible,  and  then  pol¬ 
ished.  Obtaining  a  cross-section  that  includes  the  larval  shell 
after  polishing  is  nearly  impossible  in  our  experience. 

Shells,  statoliths,  and  otoliths  are  three-dimensional  struc¬ 
tures,  and  the  materials  underlying  the  carlxmate  of  interest 
may  confound  elemental  analysis.  Based  on  our  results,  laser 
ablation  ICP-MS  analysis  may  result  in  contamination  of  natal 
habitat  signatures  by  underlying  carbonate  material  formed 
later  in  the  organism's  life.  Laser  burn-through  issues  are  espe¬ 
cially  important  in  connectivity  studies  where  it  is  often  nec¬ 
essary  to  target  an  extremely  small  area  of  a  shell  that  repre¬ 
sents  the  material  deposited  before  potential  dispersal  away 
from  natal  locations.  We  emphasize  that  any  system  should  be 
closely  scrutinized  to  assure  that  analyses  are  not  compro¬ 
mised  by  laser  burn-through. 
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Abstract 

The  chemical  composition  of  bivalve  shells  can  reflect  that  of  their  environment,  making  them  useful  indicators  of  climate,  pollution,  and 
ecosystem  changes.  However,  biological  factors  can  also  influence  chemical  properties  of  biogenic  carbonate.  Understanding  how  these  factors 
affect  chemical  incorporation  is  essential  for  studies  that  use  elemental  chemistry  of  carbonates  as  indicators  of  environmental  parameters.  This 
study  examined  the  effects  of  bivalve  shell  growth  rate  and  age  on  the  incorporation  of  elements  into  juvenile  softshell  clams,  Mya  arenaria. 
Although  previous  studies  have  explored  the  effects  of  these  two  biological  factors,  reports  have  differed  depending  on  species  and  environmental 
conditions.  In  addition,  none  of  the  previous  studies  have  examined  growth  rate  and  age  in  the  same  species  and  within  the  same  study.  We  reared 
clams  in  controlled  laboratory  conditions  and  used  solution-based  inductively  coupled  plasma  mass  spectrometry  (ICP-MS)  analysis  to  explore 
whether  growth  rate  affects  elemental  incorporation  into  shell.  Growth  rate  was  negatively  correlated  with  Mg,  Mn,  and  Ba  shell  concentration, 
possibly  due  to  increased  discrimination  ability  with  size.  The  relationship  between  growth  rate  and  Pb  and  Sr  was  unresolved.  To  determine  age 
effects  on  incorporation,  we  used  laser  ablation  ICP-MS  to  measure  changes  in  chemical  composition  across  shells  of  individual  clams.  Age 
affected  incorporation  of  Mn,  Sr,  and  Ba  within  the  juvenile  shell,  primarily  due  to  significantly  different  elemental  composition  of  early  shell 
material  compared  to  shell  accreted  later  in  life.  Variability  in  shell  composition  increased  closer  to  the  umbo  (hinge),  which  may  be  the  result  of 
methodology  or  may  indicate  an  increased  ability  with  age  to  discriminate  against  ions  that  are  not  calcium  or  carbonate.  The  effects  of  age  and 
growth  rate  on  elemental  incorporation  have  the  potential  to  bias  data  interpretation  and  should  be  considered  in  any  biogeochemical  study  that 
uses  bivalves  as  environmental  indicators. 

Published  by  Elsevier  B.V. 
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I.  Introduction 

Marine  bivalves  have  been  used  successfully  as  indicators  of 
environmental  properties  for  several  decades  (e.g.  Dodd,  1965; 
Klein  et  al..  1996;  Boisson  et  al.,  1998).  Bivalves  are  ideally 
suited  for  this  purpose  given  their  sedentary  nature  after  re¬ 
cruiting  to  the  benthos,  high  abundances,  relatively  large  sizes, 
longevity,  and  hardiness  (Phillips,  1977).  In  particular,  bivalve 
shells  have  proven  useful  for  environmental  reconstructions 
since  shell  and  ambient  water  elemental  concentrations  exhibit  a 
monotonically  increasing  relationship  (Wilbur.  1972).  This 
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relationship  is  sometimes  affected  by  environmental  properties 
such  as  temperature  and  salinity;  researchers  exploit  this 
characteristic  to  explore  environmental  conditions  that  occurred 
during  shell  development  (e.g.  Rucker  and  Valentine,  1961; 
Dodd,  1965;  Lerman,  1965;  Dodd  and  Crisp,  1982;  Pitts  and 
Wallace,  1994;  Lazareth  et  al.,  2003).  In  studies  where  bivalves 
are  used  as  indicators  of  environment,  whole  shells  or  portions 
of  individual  shells  are  analyzed  for  their  elemental  composi¬ 
tion,  and  then  related  to  spatial  or  temporal  variation  in  ele¬ 
mental  concentrations,  temperature,  or  salinity  of  the  seawater 
of  formation. 

Although  the  relationship  between  shell  composition  and 
water  chemistry  has  been  studied  in  the  past,  there  are  no 
consistent  results  across  different  species.  For  instance,  Sr:Ca  in 
molluscan  shell  has  been  reported  as  correlating  both  positively 
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(Dodd,  1965;  Smith  et  al.,  1979;  Steelier  et  al.,  1996)  and 
negatively  (Zacherl  et  al.,  2003)  to  temperature,  depending  on 
the  species  examined.  These  differences  likely  originate  in  the 
physiological  underpinnings  of  shell  formation  and  deposition, 
such  as  regulation  of  extrapallial  fluid  composition  by  the 
mantle,  which  differs  depending  on  species  (Wilbur,  1972). 
Such  contradictory  results  highlight  the  need  for  understanding 
biological  factors  that  potentially  influence  elemental  composi¬ 
tion  of  shell. 

Despite  some  uncertainty  about  the  physiological  processes 
that  drive  the  monotonic  relationship  between  shell  chemistry 
and  water  composition,  many  fields  of  study  take  advantage 
of  the  relationship  to  answer  questions  about  the  environment 
and  the  ecosystem  as  a  whole.  For  instance,  environmental 
managers  use  bivalves  as  bioindicators  of  metal  contamination 
in  coastal  habitats  by  measuring  the  trace  metal  content  of  shells 
(Lindh  et  al.,  I98X;  Bourgoin,  1990;  Fuge  et  al.,  1993;  Pearce 
and  Mann,  2006).  Another  use  of  bivalve  shell  composition  is  as 
a  recorder  of  habitat  use.  For  instance,  population  biologists  are 
exploring  the  use  of  geographically  variable  chemical  signa¬ 
tures  recorded  in  larval  shell  as  natural  tags  to  track  dispersal 
(Zacherl,  2005;  Becker  et  al.,  2007).  Shell  chemistry  also  can  be 
used  as  a  proxy  for  changes  in  water  properties.  Biological 
monitors  of  water  chemistry,  such  as  bivalves,  are  especially 
important  in  habitats  that  are  difficult  to  extensively  sample  on 
relevant  time  scales,  such  as  deep-sea  hydrothermal  vents  (Hart 
and  Blusztajn,  1998),  or  in  paleoclimatic  studies  of  salinity  or 
temperature  change  on  millennial  time  scales  (e.g.  Bourgoin 
and  Risk,  19X7). 

In  this  study,  we  were  most  concerned  with  biological  factors 
that  might  affect  elemental  incorporation  into  shell,  specifically 
the  effect  of  growth  rate  and  age.  Declines  in  elemental 
incorporation  into  the  shell’s  crystal  latice  have  been  reported  to 
occur  with  age  for  Pb  in  the  abalone  Haliotis  spp.  (Hirao  et  al.. 

1 994;  Arai  et  al..  2003)  and  for  Sr  and  Mg  in  the  bivalve  Mytilus 
spp.  (Dodd,  1965).  In  addition  to  age,  growth  rate  may  impact 
elemental  incorporation  into  shell.  Individuals  with  rapid 
growth  rates  have  been  shown  to  incorporate  higher  amounts 
of  elements  into  carbonate  compared  to  slower-growing 
individuals  of  the  same  species  (Wilbur,  1972).  This  effect 
has  been  observed  in  coral  skeletons  (Marshall  and  McC  ulloch, 
2002;  Mitsuguchi  et  al..  2003),  fish  otoliths  (Hamer  and 
Jenkins,  2007),  and  adult  bivalve  shells  (Stecher  et  al.,  1996; 
Ciillikin  et  al.,  2005;  Carre  et  al.,  2006).  These  two  biological 
factors  have  the  potential  to  decouple  relationships  between 
ambient  water  properties  and  biogenic  carbonate  composition. 
As  such,  they  must  be  understood  and  taken  into  account  during 
any  attempts  to  reconstruct  environmental  conditions  or  habitat 
use  based  on  shell  chemistry. 

Although  the  studies  above  have  explored  the  effects  of 
growth  rate  and  age  on  elemental  incorporation  into  carbonate, 
none  have  explored  the  effects  of  these  two  biological  factors  in 
the  same  species  and  within  the  same  study.  The  goal  of  this 
study  was  to  understand  the  relationships  between  growth  rate 
and  age  on  shell  chemistry  in  the  commercially  important 
softshell  clam,  Mva  arenaria.  First,  we  examined  the  effects  of 
growth  rate  on  shell  elemental  composition  by  comparing  shells 


of  clams  from  the  same  cohort,  reared  in  the  same  conditions, 
but  with  different  final  sizes  and  therefore  different  growth 
rates.  Second,  we  compared  different  regions  of  shell  within  the 
juvenile  stage  to  explore  the  effects  of  age  on  incorporation. 
Previous  studies  have  focused  on  one  or  two  elements  for 
exploring  the  relationship  between  biological  factors  and 
elemental  incorporation  into  shell  (e.g.  Hirao  et  al..  1994; 
Ciillikin  et  al.,  2005;  Carre  et  al..  2006).  Here  we  look  at  five 
elements  found  to  be  useful  indicators  of  environmental 
properties  in  previous  studies.  In  addition,  we  used  controlled 
laboratory  conditions  to  distinguish  between  biological  factors 
and  environmental  factors  that  vary  in  natural  settings,  such  as 
temperature  and  salinity.  Our  results  have  implications  for 
biogeochemical  studies  that  use  bivalve  shell  elemental  com¬ 
position  as  an  indicator  of  environmental  parameters.  If  age  or 
growth  rate  affects  incorporation  of  elements  into  carbonate, 
investigators  are  obliged  to  take  these  factors  into  consideration 
when  interpreting  environmental  variables  based  on  shell 
chemistry. 

2.  Methods 

2.1.  Clam  rearing 

Adult  M.  arenaria  with  ripe  gonads  were  obtained  from 
Cotuit,  Massachusetts  in  April  2006  and  transported  to  the 
Environmental  Systems  Laboratory  (ESL)  at  Woods  Hole 
Oceanographic  Institution,  where  they  were  placed  in  mesh 
bags  and  suspended  in  a  750  L  tank  with  filtered  seawater. 
Spawning  activity  commenced  approximately  I  h  later,  and  the 
tank  was  left  undisturbed  to  allow  spawning  to  complete  and  for 
fertilization  to  take  place.  After  5  h  the  adult  clams  were 
removed,  and  the  tank  contents  were  filtered  through  a  35  pm 
synthetic  nylon  mesh  sieve  to  concentrate  the  larvae  into  a  small 
volume  (~20  L).  We  counted  a  subsample  of  larvae  and 
obtained  a  total  estimate  of  16  million  trochophore  larvae. 

Trochophore  larvae  were  placed  into  12  L  high-density 
polyethylene  tanks  (three  tanks  per  experiment)  at  a  density  of 
approximately  45  larvae  mL  '.  Clams  whose  shells  were 
intended  for  the  growth  rate  experiment  (see  explanation  of 
experiments  below)  were  reared  in  water  with  salinity  ~22.5%<» 
at  20  °C.  Clams  whose  shells  were  intended  for  the  age 
experiment  were  reared  in  undiluted  seawater  (~30%o)  at  24  °C. 
Temperatures  and  salinities  were  within  the  range  experienced 
by  this  species  in  temperate  tidal  estuaries  and  were  chosen 
because  the  tanks  were  in  use  as  part  of  a  larger  set  of 
experiments;  no  comparisons  were  made  between  shells  of 
clams  reared  in  different  salinities.  Tanks  were  placed  into  large 
water  baths  to  maintain  experimental  temperatures.  The  20  °C 
water  bath  received  a  continuous  supply  of  20  °C  water,  which 
was  regulated  for  the  facility's  seawater  supply  line  by  large- 
volume  chillers  and  heaters.  The  water  bath  for  24  °C  tanks  was 
maintained  using  a  120-volt  tank  heater  regulated  by  a 
thermostat  (Process  Technology  Easy  Plug™  Heater  with 
Digital  Controller,  1X00  W).  Seawater  for  tanks  was  obtained 
from  the  in-house  supply  line,  which  pumps  water  from 
Vineyard  Sound  100  m  offshore  at  a  water  depth  of  4  m.  All 
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seawater  was  filtered  to  remove  particles  >1  pm  before  use. 
Salinity  was  reduced  to  22.5%o  by  adding  ultrapure  H20  to 
Vineyard  Sound  seawater.  Tank  temperatures  and  salinities 
were  measured  every  other  day,  and  water  samples  were  taken 
weekly  to  determine  ambient  elemental  ratios  over  the  course  of 
the  experiment  (see  Section  2.2). 

Clams  were  raised  for  60  days  under  treatment  conditions, 
with  complete  water  changes  every  2  days.  New  tank  water  was 
adjusted  to  the  appropriate  temperature  and  salinity  before 
clams  were  added.  Larvae  metamorphosed  into  juveniles  within 
2  weeks  of  spawning.  Larvae  and  juveniles  were  fed  a  mix  of 
live  Isochrysis  sp.  and  concentrated  algae  (Instant  Algae© 
Shellfish  Diet  1800)  one  to  two  times  daily  based  on  clearance 
rates.  All  tanks  received  the  same  food  mixture  over  the  course 
of  the  experiment.  After  60  days,  clams  ranged  in  size  from 
0.8  mm  to  7.1  mm;  they  were  removed  from  their  tanks  and 
frozen  until  cleaned  and  prepared  for  analysis.  There  were  >200 
surviving  clams  in  each  of  the  six  tanks,  and  we  analyzed  25, 
27,  and  25  individuals  from  the  three  growth  rate  experiment 
tanks  and  7,  5,  and  5  individuals  from  the  three  age  experiment 
tanks.  These  sample  sizes  were  chosen  based  on  statistical  tests 
that  determined  ideal  sizes  of  25  and  5  for  solution  and  laser 
ablation  ICP-MS,  respectively.  We  prepared  two  additional 
samples  from  each  tank  when  possible  to  account  for  losses 
during  the  shell  cleaning  process.  If  samples  remained  intact 
after  cleaning,  we  kept  the  resulting  larger  sample  size  to 
improve  statistical  power. 

2.2.  Seawater  analysis 

We  quantified  variability  in  seawater  composition  by 
measuring  elemental  ratios  in  each  tank  for  each  week  of  the 
experiment  (w  =  9  per  tank).  Thus  it  was  possible  to  distinguish 
differences  in  elemental  incorporation  due  to  variable  seawater 
chemistry  versus  ontogeny  or  growth  rate.  Samples  were 
vacuum  filtered  using  acid-washed  plastic  funnels  with  0.2  pm 
cellulose  nitrate  membrane  filters.  Samples  were  then  trans¬ 
ferred  to  acid-washed  HDPE  bottles,  acidified  to  pH  ~2  using 
ultrapure  HNO3,  and  refrigerated  for  up  to  two  months  before 
analysis.  Seawater  samples  were  prepared  for  analysis  by 
diluting  50-fold  with  ultrapure  2%  HNO3. 

Water  samples  were  analyzed  using  a  Thermo-Finnigan 
MAT  Element2  magnetic  sector  field  inductively  coupled 
plasma  mass  spectrometer  (ICP-MS).  To  correct  for  mass  bias 
and  instrument  drift,  a  2%  HNO3  blank  solution  and  CASS-4 
Nearshore  Seawater  Reference  Material  (National  Research 
Council  Canada  Certified  Reference  Material)  were  run  periodi¬ 
cally.  During  analyses  we  monitored  25Mg,  4xCa,  55Mn,  K8Sr, 
nxBa,  and  208 Pb  in  low  resolution  mode.  These  elements  were 
chosen  because  of  their  previous  use  as  environmental  tags 
(Kalish,  1989;  Campana,  1999;  Elsdon  and  Gillanders,  2003; 
Zacherl  et  al.,  2003).  Molar  ratios  of  each  element  to  Ca  were 
calculated  using  mass  bias  corrections  calculated  from  calibration 
standards.  Limits  of  detection  (LOD,  3er  of  blank)  were  as 
follows:  27  pg  g-  1  for  Mg;  7.7  pg  g-  1  for  Ca;  5.6  pg  g-  1  for  Mn; 
5.9  pg  g- 1  for  Sr;  0.42  pg  g- 1  for  Ba;  and  0.50  pg  g~  1  for  Pb.  We 
tested  for  differences  in  seawater  chemistry  among  tanks  and  over 


time  within  a  tank  using  a  two-way  analysis  of  variance 
(ANOVA).  Significant  differences  in  seawater  composition 
were  accounted  for  by  using  partition  coefficients  to  report  shell 
elemental  incorporation.  Partition  coefficients,  or  the  ratio 
between  Element:Ca  values  in  carbonate  and  ambient  water,  are 
used  in  biogeochemistry  literature  to  relate  seawater  elemental 
ratios  to  those  of  carbonate  (Lea  and  Spero,  1 992).  We  conducted 
all  analyses  for  growth  rate  and  age  effects  using  both  elemental 
ratios  and  partition  coefficients.  Both  methods  yielded  the  same 
results,  and  we  report  only  elemental  ratios  hereafter. 

2.3.  Growth  rate  experiment 

Since  all  clams  were  from  the  same  spawning  event  and  were 
allowed  to  grow  for  60  days  after  fertilization,  differences  in 
final  shell  size  were  attributed  to  different  growth  rates.  We 
recorded  the  final  length  along  the  longest  axis  of  the  shell  to  the 
nearest  0.01  mm  of  each  individual.  This  length  was  converted 
to  growth  rate  (pm  day-1)  by  assuming  a  linear  relationship 
between  size  and  age.  Typically,  the  nonlinear  Von  Bertalanfiy 
equation  is  used  to  describe  bivalve  growth  rate  (Brousseau, 
1979;  Appeldoom,  1983).  This  equation,  however,  is  most 
applicable  to  clams  much  older  than  60  days.  The  portion  of  the 
Von  Bertalanfiy  curve  which  accounts  for  clam  growth  during 
the  first  60  days  is  approximately  linear  (Von  Berlalanffy, 
1938),  therefore  the  assumption  of  linear  growth  was  appro¬ 
priate  for  our  purposes. 

Shells  were  cleaned  thoroughly  using  techniques  developed 
for  foraminiferan  tests  (Boyle,  1981)  with  modifications 
specifically  for  M.  arenaria.  The  most  notable  modification 
was  removal  of  the  reductive  cleaning  step  because  it  dissolved 
the  proteinaceous  structure  of  the  shell  so  that  it  could  not  be 
prepared  for  laser  ablation.  Clam  shells  were  placed  in 
individual  acid-washed  vials  using  acid-washed  plastic  forceps, 
and  sonicated  briefly  to  remove  organic  matter.  Individuals 
were  rinsed  three  times  with  ultrapure  water,  and  then  soaked 
for  10  min  at  80  °C  in  1%  H202  solution  buffered  in  1  N 
ultrapure  NaOH  to  remove  organic  material.  Afterward,  shells 
were  rinsed  three  times  with  ultrapure  water,  transferred  to 
clean,  acid-washed  vials,  rinsed  four  times  with  ultrapure  water, 
then  left  to  dry  overnight  in  a  laminar  flow  hood.  After  shells 
were  dry,  they  were  weighed  to  the  nearest  0. 1  mg  and  dissolved 
in  2%  ultrapure  HN03  to  achieve  a  20,000-fold  dilution  of 
calcium  carbonate  for  each  clam  based  on  weight.  All  shell 
cleaning  and  preparation  was  done  in  a  Class  100  clean  room. 

Shell  material  was  analyzed  using  ICP-MS  with  corrections 
for  mass  bias  and  instrument  drift  as  for  seawater  analyses  using 
two  solution-based  standards,  an  aragonitic  otolith  reference 
material  (Yoshinaga  et  al.,  2000)  and  the  certified  reference 
material  FEBS-1  (Sturgeon  et  al.,  2005).  The  same  elements 
noted  above  were  analyzed  in  low  resolution  mode  and 
converted  into  molar  ratios  relative  to  Ca  using  mass  bias 
corrections  calculated  from  calibration  standards.  Limits  of 
detection  were  as  follows:  0.218  pg  g-  1  for  Mg;  0.819  pg  g-  1 
for  Ca;  1 .36  pg  g- 1  for  Mn;  0.632  pg  g- 1  for  Sr;  0. 1 50  pg  g-  1 
for  Ba;  and  0.122  pg  g- 1  for  Pb.  Elemental  ratios  that  were 
more  than  two  standard  deviations  from  the  mean  of  all  sample 
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Fig.  1 .  Diagram  of  juvenile  shell  showing  positions  of  the  laser  tracks  for  the  age 
experiment.  The  distances  a  =  700  pm  and  b  =  360  pm.  Lines  ablated  along 
growth  rings  (dashed)  are  1400  pm  each. 

intensities  for  a  particular  element  (i.e.  the  global  mean)  were 
eliminated  from  further  analysis,  and  the  entire  shell  was 
eliminated  if  three  or  more  elemental  ratios  met  the  exclusion 
criteria.  We  removed  17/380  (4.5%)  of  the  elemental  ratios  for 
the  age  experiment,  plus  all  data  from  one  entire  clam  shell, 
before  further  analyses.  There  was  no  evidence  of  matrix  effects 
for  any  element  based  on  plots  of  E:Ca  versus  4XCa  intensity. 

To  determine  whether  growth  rate  affected  elemental 
incorporation  into  shell,  ideally  we  would  pool  data  from 
each  replicate  tank  of  the  growth  rate  experiment  to  increase 
statistical  power.  However  because  clam  sizes  differed  among 
tanks,  we  first  performed  regression  analysis  on  growth  rate 
versus  elemental  ratios  for  each  tank  separately  (analysis  of 
covariance).  We  then  compared  slopes  among  tanks  to 
determine  the  validity  of  pooling  using  equations  from  Zar 
(1999).  In  the  case  where  slopes  were  not  significantly  different 
(i.e.  F<Fcn{),  we  calculated  the  regression  slope  for  all  three 
tanks  combined. 

Our  method  assumes  a  linear  relationship  between  growth 
rate  and  elemental  incorporation  into  shell.  If  the  relationship  is 
instead  nonlinear,  then  slope  differences  among  tanks  may 
result  from  the  different  ranges  in  growth  rates  represented  in 
each  tank.  To  test  for  this  possibility,  we  repeated  the  regression 
analysis  described  above  with  the  data  selected  to  include  only 
the  range  in  growth  rates  observed  in  all  three  tanks  (0.040- 
0.064  mm  day-1).  If  a  nonlinear  relationship  was  responsible 


for  slope  differences  among  tanks,  we  would  expect  those 
differences  to  disappear  when  only  selected  data  are  used. 

2.4.  Age  experiment 

All  clams  were  of  the  same  age  (60  days)  when  the  experiment 
was  terminated.  To  determine  variability  in  elemental  incorpora¬ 
tion  with  age,  we  measured  and  compared  elemental  composition 
at  different  growth  rings  along  the  shell  of  an  individual  clam, 
which  correspond  to  different  clam  ages.  Growth  rings  located 
closer  to  the  umbo  were  laid  earlier  in  the  clam’s  life,  and  rings 
further  from  the  umbo  were  laid  more  recently.  We  chose 
individuals  that  were  of  similar  size  to  minimize  the  effects  of 
differences  due  to  growth  rate.  Exploratory  statistics  showed  that 
there  were  no  significant  correlations  between  elemental 
composition  and  growth  rate  for  any  of  the  individuals  analyzed 
for  the  age  experiment,  allowing  us  to  interpret  the  effects  of  age 
without  the  confounding  effects  of  growth  rate. 

Shells  were  cleaned  as  in  the  growth  rate  experiment,  and 
their  lengths  were  recorded  to  the  nearest  0.01  mm.  Individual 
shell  valves  were  then  mounted  on  glass  slides  using  Devcon© 
Super  Glue.  We  were  not  able  to  polish  shells  to  an  even  plane 
for  analysis  as  this  would  have  resulted  in  excessive  sample  loss 
due  to  the  combined  effects  of  the  shells’  curvature  and  thinness 
(<500  pm).  Shell  material  was  analyzed  using  ICP-MS  coupled 
to  a  New  Wave  Research  UP213  laser.  We  ablated  material  by 
tracking  the  laser  along  five  concentric  growth  increments  of 
juvenile  clam  shells  radiating  out  from  the  umbo.  Each  ring 
ablated  was  1400  pm  in  length  and  traced  one  of  the  growth 
increments  of  the  shell.  Measurements  were  taken  at  the  same 
distances  from  the  umbo  for  all  shells,  with  the  assumption  that 
the  resulting  measurements  represented  similar  ages  among 
individuals.  This  assumption  was  valid  since  growth  rates  did 
not  differ  significantly.  The  first  ring  measured  was  700  pm 
away  from  the  umbo  as  measured  along  the  axis  of  growth,  and 
the  remaining  four  rings  measured  were  spaced  360  pm  from 
one  another  (Fig.  1 ).  The  laser  was  set  to  80%  output  (0. 1 2  mJ 
per  pulse),  with  a  10  Hz  repetition  rate,  30  pm  spot  size,  and  a 
scan  speed  of  1 0  pm  s~ 1 . 

Vaporized  material  from  the  ablation  was  transported  via  a 
helium  gas  stream  to  the  dual-inlet  quartz  spray  chamber  where 
it  was  mixed  with  1%  HNO^  aerosol  from  a  self-aspirating  PFA 
20  pL  min  1  nebulizer.  The  analyte  was  then  transported  to  the 


Table  1 

Mean  temperature,  salinity,  and  dissolved  ambient  seawater  elemental  ratios  to  calcium  (±SE)  for  tanks  from  both  experiments 


Temp 

Salinity 

Mg:Ca 

Mn:Ca 

Sr:Ca 

Ba:Ca 

Pb:('a 

°C 

%o 

(mol  moF ') 

(pmol  mol" ') 

(mmol  mol"  ) 

(pmol  mol" 1 ) 

(nmol  mol"1) 

Growth  rate  experiment 

Tank  1  19.4  ±0.04 

Tank  2  I9.4±0.04 

Tank  3  19.5  ±0.04 

22.4±0.07 

22.4  ±0.05 

22.3  ±0.05 

5. 09  ±0.04 

5.09  ±0.03 

5.05  ±0.02 

1 8.6±4.94 

1 3.6±  3.55 

1 4.7  ±  3.98 

8.56  ±0.03 

8.64±0.03 

8.57±0.03 

1 4.9  ±  1.38 

15.1  ±  1.41 

15.1  ±  1.19 

15.3  ±4. 12 

9.37±2.27 

9.32  ±  1.53 

Age  experiment 
Tank  1 

Tank  2 

Tank  3 

23.9±0.05 

23.8±0.05 

23.9±0.05 

29.5  ±0.02 

29.7  ±0.06 

29.5  ±0.25 

5.03  ±0.01 

5.03  ±0.02 

5.01  ±0.03 

1 2.8  ±  3.26 

1 1  6±3.66 

1 2.6  ±  2.65 

8.58  ±0.03 

8.61  ±0.03 

8.61  ±0.04 

14.9±0.68 

1 4.5  ±  0. 1 7 

15. 6  ±1.03 

11.0  ±2.04 

8.31  ±1.91 
9.45  ±2.26 
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Table  2 


Results  of  ANOVA  testing  for  differences  of  seawater  elemental  ratios  due  to  week  of  experiment  or  tank 


Source  DF 

Mg:Ca 

Mn:Ca 

Sr:Ca 

Ba:Ca 

Pb:Ca 

MS 

F 

P 

MS 

F 

P 

MS 

F 

P 

MS 

F 

P 

MS 

F 

P 

Growth  rate  experiment 

Week  8 

0.017 

9.441 

0.000 

479.173 

18.409 

0.000 

0.012 

3.303 

0.020 

17.085 

1.405 

0.267 

186.687 

6.899 

0.001 

Tank  2 

0.003 

1.769 

0.202 

2.699 

0.104 

0.902 

0.011 

3.101 

0.073 

0.211 

0.017 

0.983 

61.525 

2.286 

0.134 

Error  1 6 

0.002 

26.029 

0.004 

12.164 

26.917 

Age  experiment 

Week  8 

0.006 

1.884 

0.134 

183.184 

5.282 

0.002 

0.011 

1.625 

0.194 

4.193 

0.878 

0.555 

82.594 

4.936 

0.003 

Tank  2 

0.001 

0.485 

0.624 

3.680 

0.106 

0.900 

0.001 

0.174 

0.842 

2.713 

0.568 

0.578 

16.065 

0.960 

0.404 

Error  16 

0.003 

34.681 

0.007 

4.778 

16.732 

Bold  F  statistics  and  p  values  are  significant. 


ICP-MS  via  an  argon  carrier  gas.  We  corrected  for  mass  bias  and 
instrument  drift  using  standards  as  in  the  growth  rate 
experiment.  Elements  were  measured  in  medium  resolution 
mode,  and  molar  ratios  of  each  element  to  4xCa  were  calculated 
using  mass  bias  corrections  as  above.  LOD  were  as  follows: 
0. 1 53  pg  g-  1  for  Mg;  0.58 1  pg  g~  1  for  Ca;  0.836  pg  g”  1  for 
Mn;  1 .04  pg  g“  1  for  Sr;  0.095  pg  g~ 1  for  Ba;  and  0. 1 69  pg  g~  1 
for  Pb.  Elemental  ratios  that  were  more  than  two  standard 
deviations  from  the  global  mean  were  eliminated  from  further 
analysis,  and  results  for  an  entire  growth  ring  were  removed  if 
three  or  more  of  the  elemental  ratios  met  the  exclusion  criteria. 
We  removed  14/450  (3.1%)  of  the  elemental  ratios  for  the  age 
experiment,  plus  all  data  from  two  entire  clam  shells,  before 
analyses.  There  was  no  evidence  of  matrix  effects  for  any 
element  based  on  plots  of  E/Ca  versus  4xCa  intensity.  To  test  for 
differences  in  elemental  incorporation  due  to  age,  we  compared 
shell  composition  between  growth  rings  within  a  clam  using 
multivariate  repeated  measures  ANOVA  (Winer,  1971;  Barci- 
kowski  and  Robey,  1984).  Data  were  log-transformed  to 


achieve  a  normal  distribution;  log-transformed  data  were  used 
for  the  plots  and  for  all  regression  calculations. 

3.  Results 

3.1.  Seawater  analyses 

Although  seawater  for  all  tanks  originated  from  the  same 
source,  ambient  elemental  composition  varied  over  time  for 
tanks  in  both  experiments  (Tables  1  and  2).  Temporal  variability 
in  seawater  chemistry  for  the  growth  rate  experiment  was  not 
expected  to  affect  our  analyses  since  we  compared  only  the  bulk 
shell  composition  with  the  time-averaged  water  composition. 
For  the  age  experiment,  however,  temporal  variability  might 
affect  our  analyses  since  we  measured  areas  of  shell  that 
correspond  to  specific  time  periods.  Results  from  ANOVA  and 
plots  of  tank  elemental  ratios  over  the  course  of  the  age 
experiment  showed  that  tanks  differed  over  time  for  Mn:Ca  and 
Pb:Ca  (Table  2,  Fig.  2).  Elemental  ratios  for  weeks  13  were 


A 


Week 


B  C 


Week 


Fig.  2.  Seawater  elemental  ratios  over  9  weeks  for  each  tank  in  the  age  experiment  for  (A)  Mg:Ca,  (B)  Mn:Ca,  (C)  Sr:Ca,  (D)  Ba:Ca,  and  (E)  Pb:Ca  (tank  1  =  circles; 
tank  2  =  squares;  tank  3  =  crosses). 
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Fig.  3.  Elemental  ratios  of  individual  clam  shells  for  Mya  arenaria  growing  at  different  rates  for  (A)  Mg.Ca,  (B)  Mn:Ca,  (C)  Sr:Ca,  (D)  Ba:Ca.  and  (E)  Pb:C'a  Data  are 
from  three  replicate  tanks  (tank  1  =  circles;  tank  2  =  squares;  tank  3  =  stars).  Regression  lines  are  for  each  tank  separately:  tank  1  =  solid;  tank  2  =  dotted;  tank  3 
dashed  and  dotted. 


higher  than  for  weeks  6-9  for  Mn,  and  week  3  had  higher  Pb:Ca 
than  other  weeks  (/?<0.0I,  Tukey- Kramer  pairwise  tests). 

3.2.  Growth  rate  effects 

We  used  solution-based  ICP-MS  analysis  of  juvenile  clam 
shells  from  three  replicate  tanks  (w  25,  27,  25)  to  test  the 
effects  of  growth  rate  on  elemental  incorporation  into  shell 
(Fig.  3).  Growth  rates  ranged  from  0.014  to  0.1 19  mm  day~  1 
(0.050±0.023  mm  day-1,  mean±SD).  To  determine  the  rela¬ 
tionship  between  shell  elemental  incorporation  and  growth  rate, 
we  plotted  elemental  ratios  versus  growth  rates  for  each  tank 
separately  (Fig.  3).  Values  of  R 2  ranged  from  0.005  to  0.634, 
with  7  of  the  15  regression  slopes  significantly  differing  from 
zero  (p<  0.05,  Table  3).  There  was  a  general  trend  of  decreasing 
elemental  ratios  with  increasing  growth  rate,  with  the  exception 
of  Sr:Ca  and  Pb:Ca  from  tank  3  (Fig.  3). 

Slopes  significantly  differed  among  replicate  tanks  for  two  of 
the  five  elemental  ratios  ( p < 0.00 1 ,  Table  4).  For  Sr:Ca,  the 
slope  of  tank  3  differed  from  the  other  two,  and  for  Pb:Ca  the 
slopes  of  tanks  2  and  3  differed.  As  a  result  of  the  differences 
among  tanks  for  these  two  elements,  we  were  not  able  to  draw 
conclusions  from  the  regression  statistics  for  combined  tanks.  For 
the  remaining  three  elemental  ratios  (Mg:Ca,  Mn:Ca,  Ba:Ca), 
slopes  did  not  differ  significantly  among  tanks,  so  we  calculated 
the  common  slope  for  each  element  with  equation  18.30  in  Zar 
( 1999)  using  data  from  all  three  tanks  (Table  4). 

We  performed  the  same  regression  analyses  on  selected  data 
to  test  whether  significant  differences  in  slopes  were  due  to  a 
nonlinear  relationship  between  growth  rate  and  elemental 
incorporation  into  shell.  However  data  truncation  resulted  in 
the  elimination  of -50%  of  data  points  (n=  10,  9,  9),  and  as  a 
consequence  we  lost  statistical  power  and  were  not  able  to 


detect  significant  relationships  for  any  tanks  or  elements. 
We  could  not,  therefore,  use  the  selected  data  to  pool  values 
from  replicate  tanks,  and  the  issue  of  nonlinearity  remains 
unresolved. 

3.3.  Age  effects 

We  analyzed  shell  along  five  growth  rings  of  juvenile  clams 
from  three  replicate  tanks  (/?  =  7,  5,  5)  to  determine  whether 
elemental  incorporation  varied  as  the  individual  aged  (Fig.  4). 

Table  3 


Results  of  regression  analysis  for  growth  rate  experiment  for  individual  tanks 


Source 

ft. 

A 

f? 

F 

P 

Mg:Ca 

Tank  1 

-6.537 

-0.214 

0.089 

2.057 

0.166 

Tank  2 

-5  696 

-0.431 

0.434 

18.374 

0.0003 

Tank  3 

-6.896 

-0.150 

0.070 

1.662 

0.211 

Mn:Ca 

Tank  1 

-11.311 

0.084 

0.005 

0.103 

0.751 

Tank  2 

-  10.510 

-0.079 

0.004 

0.106 

0.766 

Tank  3 

-  10.124 

-0.338 

0.109 

2.682 

0.116 

Sr:Ca 

Tank  1 

6.098 

-0.031 

0.120 

3.010 

0.097 

Tank  2 

-5.775 

-0.108 

0.518 

22.558 

0.0001 

Tank  3 

-6.410 

0.0583 

0.199 

5.483 

0.029 

Ba:Ca 

Tank  1 

-12.245 

-0.303 

0.243 

7.051 

0.014 

Tank  2 

- 1 1 .629 

-0.469 

0.606 

33.861 

0.000 

Tank  3 

- 12.710 

-0.255 

0.566 

27.619 

0.000 

Pb:Ca 

Tank  1 

-17.262 

-0.602 

0.144 

3.363 

0.082 

Tank  2 

-13.570 

-1.608 

0.634 

41.490 

0.000 

Tank  3 

-21.435 

0.387 

0.032 

0.737 

0.400 

The  coefficients  ft,  and  ft  correspond  to  the  regression  equation,  where y  =ft>+ft.v. 
Bold  F  statistics  and  p  values  are  significant. 
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Table  4 

Results  of  statistic  testing  for  differences  among  regression  functions  of  replicate 
tanks 


Elemental  ratio 

F 

P 

Pi 

Mg:Ca 

1.320 

0.273 

-0.271 

Mn:Ca 

0.745 

0.465 

0.078 

Sr:Ca 

12.936 

<0.0005 

-0.032 

Ba:Ca 

1.439 

0.244 

-0.347 

Pb:Ca 

8.010 

<0.001 

-0.692 

Bold  F  statistics  are  significant  and  indicate  that  slopes  are  different  among 
tanks,  making  it  invalid  to  pool  data. 


Repeated  measures  multivariate  ANOVA  indicated  that  there 
were  significant  differences  in  elemental  ratios  among  growth 
rings  for  Mn:Ca,  Sr:Ca  and  Ba:Ca  (Table  5).  In  all  three  cases, 
the  first  growth  ring  (0.70  mm  away  from  the  umbo)  differed 
from  the  remaining  four  growth  rings  for  individuals  in  one  or 
more  replicate  tanks  (Tukey- Kramer  pairwise  comparisons). 
There  were  also  significant  differences  among  measurements  of 
clam  shells  from  the  same  tank  for  Mn:Ca,  Ba:Ca,  and  Pb:Ca. 
Finally,  we  detected  a  tank  effect  for  one  elemental  ratio;  mean 
Ba:Ca  of  shells  in  tank  1  was  higher  than  in  tanks  2  or  3. 

Although  clam  shells  of  similar  size  were  chosen  for  the  age 
experiment,  total  lengths  ranged  from  3.36  to  6.45  mm  (5.32  ± 
0.77  mm,  mean±SD),  as  measured  along  the  shell  approxi¬ 
mately  perpendicular  to  the  axis  of  growth.  The  difference  in 
size  after  60  days  resulted  in  a  range  of  growth  rates  from  0.056 
to  0. 107  mm  day-  1  (0.089±0.013  mm  day- ',  mean±SD).  We 
tested  for  a  correlation  between  growth  rate  and  elemental 
incorporation  for  each  clam  (averaged  over  growth  rings) 
analyzed  in  the  age  experiment  and  found  none  ( R 2  ranged 
from  0.0007  to  0.124;  p  values  ranged  from  0.15  to  0.91).  This 
lack  of  correlation  suggests  that  the  range  of  growth  rates  for  clam 


shells  was  sufficiently  small  so  as  not  to  impact  elemental 
incorporation. 

4.  Discussion 

4.  /.  Growth  rate  effects 

Previous  studies  have  suggested  that  the  rate  of  calcium 
carbonate  crystal  formation,  which  is  closely  tied  to  growth  rate, 
influences  elemental  incorporation  in  bivalves  (Stccher  et  al., 
1996;  Gillikin  et  al.,  2005;  Carre  et  al.,  2006).  Intuitively,  a 
higher  growth  rate  might  be  expected  to  result  in  more  crystal 
defects  during  shell  formation:  increased  active  transport  of  Ca2+ 
molecules  into  the  extrapallial  fluid  would  lead  to  higher  rates  of 
inclusion  of  non-Ca2+  ions  that  are  of  similar  size  and  charge 
(Wilbur  and  Saleuddin,  1983).  Indeed,  higher  growth  rates  have 
been  reported  as  corresponding  to  increased  inclusion  of  Mg,  Mn, 
and  Ba  in  fish  otoliths  (Bath  Martin  and  Thorrold,  2005;  Hamer 
and  Jenkins,  2007),  and  bivalve  shells  (Stecher  et  al.,  1996;  Carre 
et  al.,  2006).  However  in  our  study,  we  found  that  elemental  ratios 
of  Mg:Ca,  Mn:Ca,  and  Ba:Ca  were  all  negatively  correlated  to 
growth  rate  in  M.  arenaria.  One  explanation  for  our  disparate 
results  may  be  that  the  organism’s  physiological  ability  to 
discriminate  between  Ca2+  and  other  ions  improves  with  size. 
Clams  with  higher  growth  rates  would  have  reached  the  threshold 
size  for  improved  discrimination  for  Ca2+  sooner  than  those  with 
slower  growth  rates.  As  a  consequence,  proportionally  more  shell 
would  have  been  laid  with  lower  elemental  ratios  to  calcium  over 
our  two-month  study,  resulting  in  the  negative  correlation  that  we 
observed.  Although  such  an  age  effect  has  not  been  observed 
in  M.  arenaria ,  previous  studies  of  molluscs  have  reported 
decreased  elemental  incorporation  with  size  (Dodd,  1965;  Hirao 
et  al.,  1994;  Arai  et  al.,  2003). 


Fig.  4.  Mean  elemental  ratios  for  Mya  arenaria  from  three  replicate  tanks  for  shell  measurements  taken  at  different  distances  from  the  umbo  for  (A)  Mg:Ca,  (B)  Mn:Ca, 
(C)  Sr:Ca,  (D)  Ba:Ca,  and  (E)  Pb:Ca  (tank  1  =  circles;  tank  2  =  squares;  tank  3  =  stars).  Shaded  regions  are  ±  1  SD  of  mean. 
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Table  5 

Results  of  multivariate  repeated  measures  ANOVA  testing  for  differences  between  growth  rings  in  three  replicate  tanks 


Mg:Ca 

Mn:Ca 

Sr:Ca 

Ba:Ca 

Pb:Ca 

Univariate 

DF 

MS 

F 

P 

MS 

F 

P 

MS 

F 

P 

MS 

F 

P 

MS 

F 

P 

Among  tanks 
Within  tanks 

2 

0.33 

0.93 

0.42 

0.01 

0.68 

0.53 

0.07 

1.92 

0.19 

9.11 

25.29 

0.0001 

10.1 

2.36 

0.13 

Among  rings 

2 

0.22 

0.86 

0.49 

0.01 

3.25 

0.02 

0.05 

2.65 

0.09 

0.72 

4.44 

0.01 

5.02 

3.48 

0.03 

Ring  x  Tank 

8 

0.14 

0.53 

0.79 

0.01 

1.33 

0.25 

0.03 

1.80 

0.16 

0.25 

0.53 

0.21 

0.83 

0.80 

0.74 

Multivariate 

DF 

A 

F 

P 

A 

F 

P 

A 

F 

P 

A 

F 

P 

k 

F 

P 

Among  rings 

4 

0.70 

1.06 

0.42 

0.38 

4.05 

0.03 

0.27 

5.41 

0.02 

0.31 

4.44 

0.04 

0.52 

2.56 

0.10 

Ring  x  Tank 

8 

0.74 

0.40 

0.91 

0.49 

1.07 

0.42 

0.24 

2.06 

0.10 

0.47 

0.92 

0.52 

0.49 

1.18 

0.35 

Bold  F  statistics  and  p  values  are  significant. 


Although  every  attempt  was  made  to  assure  environmental 
conditions  were  consistent  across  tanks,  we  found  a  wide  range 
of  growth  rates  for  clams  among  the  three  replicate  tanks.  Our 
analyses  indicated  no  significant  differences  in  water  chemistry, 
temperature,  or  salinity  (Tables  1  and  2),  however  there  may  be 
additional  factors  we  did  not  account  for  that  are  the  source  of 
the  variability  in  growth  rates.  For  instance,  different  biological 
conditions  were  present  in  each  of  the  tanks  owing  to  a  variety 
of  processes  potentially  occurring,  such  as  algal  growth  and 
microbial  activity.  Although  we  sampled  the  seawater  weekly  to 
quantify  the  changes  in  its  chemistry  over  time,  we  did  not 
attempt  to  identify  or  quantify  biological  activity.  This  biolo¬ 
gical  activity  might  in  turn  affect  clam  growth  and  elemental 
incorporation  into  shell.  For  example,  increases  in  food  supply 
due  to  algal  growth  would  result  in  faster  growth  rates,  or  the 
presence  of  additional  oxygen-consuming  organisms  might 
result  in  decreased  oxygen  supply  to  the  clams  and  therefore 
reduced  metabolic  and  growth  rates. 

The  most  notable  tank  effect  was  seen  in  the  relationships 
between  growth  rate  and  Sr:Ca  and  Pb:Ca;  we  found  positive  or 
negative  correlations  depending  on  tank.  Our  mixed  results  are 
particularly  interesting  for  Sr:Ca  since  previously  reported 
relationships  indicated  positive  correlations  in  bivalves  (Stecher 
et  al.,  1996;  Ciillikin  et  al.,  2005)  and  corals  (Weber,  1973).  Sr:Ca 
correlations  to  growth  rate  from  otoliths  are  mixed  (Kalish,  1 989; 
Arai  et  al.,  1 996;  Bath  et  al.,  2000;  Martin  et  al.,  2004)  but  studies 
generally  report  negative  correlations  between  Sr:Ca  and  growth 
rate  (Sadovy  and  Severin,  1992,  1994;  Hamer  and  Jenkins,  2007; 
Lin  et  al.,  2007).  Since  Sr  ions  are  of  the  same  charge  as  Ca  ions 
and  only  slightly  larger,  they  are  substituted  directly  into  the 
aragonite  crystal  lattice  (Speer,  1983).  As  a  result,  Sr:Ca  ratios 
should  decrease  with  calcification  rate  (i.e.  growth  rate)  since 
higher  Ca  concentrations  in  the  extrapallial  fluid  would  dilute  Sr 
ions  (Sinclair,  2005).  Consequently  we  expected  to  find  negative 
correlations  for  Sr:Ca  and  growth  rate  in  all  of  the  replicate  tanks, 
and  the  mixed  results  suggest  that  other  factors  are  influencing  Sr 
incorporation  into  bivalve  shell. 

We  cannot  attribute  our  results  to  variable  Sr:Ca  available  in 
seawater  since  we  found  no  significant  differences  in  Sr:Ca 
among  tanks  (Table  2).  The  more  likely  cause  is  the  different 
ranges  in  growth  rate  depending  on  tank.  Tank  3  growth  rates 
were  higher  on  average  than  those  of  tanks  I  and  2,  and  tank  3  is 
the  only  replicate  tank  with  clams  having  growth  rates  higher  than 


75  pm  day  1 .  Indeed,  shells  from  tank  3  had  Sr:Ca  and  Pb:Ca  that 
negatively  correlated  with  growth  rate,  while  shells  from  the  other 
two  tanks  positively  correlated.  There  may  be  different 
physiological  processes  operating  as  the  clam  grows  larger, 
causing  a  shift  in  the  correlation  between  growth  rate  and 
elemental  ratios  with  size.  For  instance,  Carre  et  al.  (2006) 
found  that  curved  shell  sections  in  bivalves  had  higher  Sr:Ca 
than  flat  sections.  This  occurs  because  there  is  more  organic 
matrix  in  the  curved  sections  of  shell  than  in  flat  sections,  and 
therefore  also  more  binding  sites  for  Ca2+  and  its  competing 
ions  (Rosenberg  and  Hughes.  1991).  As  the  clam  grows, 
proportionally  more  shell  is  composed  of  flat  sections,  and 
whole-shell  analysis  of  a  larger  individual  would  have  lower  Sr:Ca 
and  Pb:Ca  than  a  smaller  individual. 

In  addition  to  a  wide  range  of  growth  rates  for  clams  among 
tanks,  we  found  a  wide  range  within  tanks  as  well.  Explanations 
for  the  observed  range  within  tanks  cannot  be  attributed  to 
biological  activity  or  chemical  differences  since  these  factors 
would  affect  all  of  the  individuals  of  a  particular  tank.  The 
source  of  variability  within  tanks  is  therefore  likely  to  be  due  to 
physiology  at  the  individual  clam  level.  Metabolic  rate  and 
growth  rate  are  closely  linked,  and  differences  in  metabolic  rate 
among  clams  might  result  in  different  growth  rates  (Bayne  and 
Newell,  1983;  Rosenberg  and  Hughes,  1991).  Variable  meta¬ 
bolic  rates,  and  therefore  growth  rates,  may  originate  from 
genetic  variability.  If  the  genetic  types  of  larvae  produced  in  our 
laboratory  spawning  varied  widely  in  their  metabolic  rates,  then 
we  might  expect  to  see  differences  among  individuals  within  a 
given  tank. 

4.2.  Age  effects 

Based  on  significant  differences  among  growth  rings  within 
individual  clams,  age  significantly  influenced  the  incorporation 
of  Mn,  Sr  and  Ba  into  M.  arenaria  juvenile  shell.  This  result  is 
most  likely  attributable  to  increased  variability  in  elemental 
ratios  with  decreasing  distance  from  the  umbo.  Visual 
inspection  of  data  plots  indicates  that  variability  in  elemental 
ratios  is  highest  near  the  umbo  (Fig.  4).  Indeed,  statistical 
analyses  confirm  that  all  significant  differences  among  growth 
rings  disappear  when  the  first  growth  ring  measured  is  removed 
from  analyses.  This  result  may  be  an  artifact  of  the  laser  ablation 
method  used.  Although  laser  parameters  are  consistent  over  all 
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measurements,  the  time  period  that  ablated  shell  material 
represents  is  not  known.  For  instance,  if  a  clam  produces  more 
material  per  day  at  week  4  compared  to  previous  weeks,  then 
measurements  taken  before  week  4  will  represent  a  larger  time 
period  sampled,  and  therefore  may  result  in  larger  variability  in 
the  elemental  ratios  (see  Elsdon  and  Gillanders,  2003  for  a 
complete  discussion  of  problems  associated  with  the  analysis  of 
elements  in  calcified  structures).  This  hypothesis  contradicts  our 
assumption  of  linear  growth  over  the  course  of  our  experiment, 
however  the  data  are  not  sufficiently  conclusive  to  render  our 
assumption  false.  Another  explanation  for  increased  variability 
in  the  first  growth  ring  measured  might  be  that  seawater  varied 
more  at  the  beginning  of  the  experiment  (Fig.  2).  However  this 
hypothesis  is  not  testable  since  we  were  not  able  to  match  shell 
growth  rings  to  particular  experimental  weeks,  and  therefore  to 
specific  water  chemistry.  This  is  because  growth  patterns  in 
molluscan  shell  are  a  result  of  complex  interactions  between 
physiology  and  environment,  especially  during  early  shell 
formation,  preventing  accurate  estimates  of  the  shell  formation 
timeline  in  daily  or  weekly  increments  (see  Lutz  and  Rhodes, 
1980  for  a  complete  discussion).  Furthermore,  based  on 
observations  over  the  course  of  the  experiment,  the  first  growth 
ring  measured  was  likely  accreted  around  week  four,  at  which 
point  seawater  chemistry  variability  was  in  decline. 

Our  results  for  Mn:Ca  are  at  least  in  part  attributable  to 
changes  in  seawater  Mn:Ca  over  time.  Mn:Ca  significantly 
differed  in  tank  seawater  over  the  course  of  the  experiment 
(Table  2,  Fig.  2).  There  are  two  possible  explanations  for  this 
pattern.  First,  organic  ligands  may  be  in  higher  quantities  later  in 
the  experiment  due  to  excretion  by  clams  or  biological  activity 
relating  to  food  input;  these  ligands  may  bond  to  Mn  and  remove 
it  from  solution  (Libes,  1 992).  We  would  detect  a  drop  in  the  Mn: 
Ca  ratios  in  seawater  over  time  as  a  result,  even  though  the  input 
water  Mn:Ca  is  unchanged.  Another  explanation  is  that  some 
event  occurred  between  seawater  samples  taken  at  weeks  4  and  5 
that  caused  a  decrease  in  Mn:Ca  in  the  seawater.  The  unknown 
event  might  have  been  natural,  affecting  the  source  of  seawater 
to  the  supply  line,  or  it  may  have  occurred  within  the  supply 
infrastructure  (e.g.  the  pipes  were  flushed  to  remove  a  blockage, 
causing  seawater  chemistry  to  change).  Clams  grew  to  ~  1  mm  in 
length  (as  measured  perpendicular  to  the  axis  of  growth)  until 
around  week  4,  which  corresponds  to  approximately  0.8  mm  in 
width  along  the  axis  of  growth.  It  is  therefore  likely  that  only  the 
first  growth  ring  measured,  located  0.70  mm  from  the  umbo, 
would  have  been  affected  by  the  elevated  Mn:Ca  ratios  seen  in 
the  seawater.  Indeed,  the  highest  Mn:Ca  level  for  seawater  was 
during  week  3  in  tank  1 ,  where  clam  shells  also  had  the  highest 
two  Mn:Ca  ratios  measured,  for  growth  rings  located  0.70  and 
1 .06  mm  away  from  the  umbo. 

4.3.  Conclusions  and  future  directions 

In  this  study,  we  found  that  both  growth  rate  and  age 
significantly  affect  elemental  incorporation  into  the  shell  of 
M.  arenaria.  Growth  rate  significantly  correlated  negatively 
with  elemental  ratios  for  7  of  the  15  analyses,  which  was 
surprising  since  previously  published  studies  tended  to  report 


positive  correlations.  There  were  significant  tank  effects  in  Sr:Ca 
and  Pb:Ca  correlations  to  growth  rate,  potentially  due  to  disparate 
biological  activity  occurring  in  the  three  replicate  tanks.  We  found 
that  growth  rates  varied  widely  within  tanks,  indicating  that  there 
is  a  physiological  factor  that  we  did  not  control  for  in  our  ex¬ 
periment,  potentially  originating  at  the  genetic  level.  Our  analysis 
of  the  effects  of  age  on  elemental  incorporation  into  M.  arenaria 
shell  revealed  higher  variability  closer  to  the  umbo  than  further 
away,  which  caused  significant  differences  with  age  for  clams. 
When  the  first  growth  ring  measured  was  removed  from  analyses, 
these  significant  differences  disappeared.  We  hypothesize  that  the 
age  effect  is  primarily  a  result  of  the  shell’s  more  pronounced 
curvature  near  the  umbo,  which  affects  physiological  processes 
such  as  uptake  of  Ca2+  relative  to  other  ions  of  similar  size  and 
charge. 

The  effects  of  size  on  elemental  incorporation  have  the 
potential  to  bias  data  interpretation  in  any  biogeochemical 
study  that  uses  bivalves  as  indicators  of  environment.  Juvenile 
bivalves  are  ideal  candidates  for  environmental  indicators:  pre¬ 
dictable  spawning  behavior  allows  one  to  be  sure  what  season 
their  shell  is  laid,  and  they  are  more  susceptible  to  some  envi¬ 
ronmental  contaminants,  making  them  an  early  indicator  of 
problems  that  might  not  yet  affect  adults.  However  studies  that 
use  juveniles  in  this  capacity  must  be  careful  to  take  into  account 
the  variability  in  shell  composition  with  age  and  growth  rate. 
Variability  in  shell  chemistry  could  be  interpreted  as  indicating 
shifting  environmental  conditions  when  in  fact  composition  may 
reflect  individual  or  ontogenetic  variability  in  physiological 
processes.  As  our  data  show,  clams  reared  in  identical  environ¬ 
ments  may  incorporate  different  elemental  signatures  in  their 
shells  depending  on  their  age  or  growth  rate.  Care  must  be  taken  to 
understand  the  potentially  confounding  effects  of  physiology  for 
this  and  other  bivalve  species  before  making  inferences  about  the 
environment  based  on  shell  chemistry. 

This  study  was  the  first  to  explicitly  examine  the  effects  of 
variable  growth  rates  and  age  in  bivalves  on  incorporation  of  a 
suite  of  elements  into  shell,  and  in  the  same  species  during  the 
same  study.  Although  the  variability  we  found  is  not  cause  for 
dismissing  bivalves  as  useful  environmental  indicators,  further 
studies  should  focus  on  the  cause  of  variable  incorporation  within 
a  cohort,  especially  with  relation  to  the  variables  we  were  not  able 
to  control  for  in  this  set  of  experiments.  The  bioavailable  levels  of 
elements  in  the  seawater  should  be  measured  to  ascertain  whether 
the  source  of  variability  in  ratios  is  due  to  water  chemistry  or 
physiology.  The  sizes  of  clams  should  be  tracked  and  recorded 
over  the  entire  experiment  to  assess  individual  growth  rates  rather 
than  averages,  and  to  identify  exactly  when  individual  growth 
rings  are  accreted  so  they  can  be  matched  to  the  correct  water 
sample  taken.  Parental  analysis  would  be  useful  to  establish 
genetic  contributions  of  individuals  and  to  determine  if  certain 
genetic  combinations  lead  to  higher  or  lower  growth  rates, 
explaining  the  variability  that  we  saw  within  tanks. 
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Chapter  4 


Temperature  and  salinity  effects  on 
elemental  uptake  in  the  shells  of 
larval  and  juvenile  softshell  clams 
(My  a  arenaria) 

4.1  Introduction 

Most,  marine  benthic  invertebrate  life  cycles  include  a  planktonic  larval  phase  that 
facilitates  dispersal  among  adult  populations  (Thorson,  1950).  Connectivity,  or  the 
degree  to  which  geographically  separated  populations  exchange  individuals,  is  an  im¬ 
portant  factor  in  the  spatial  population  dynamics  of  many  marine  organisms  (Moila- 
nen  and  Nieminen,  2002).  Understanding  connectivity  in  marine  benthic  populations 
is  important  because  of  the  role  spatial  dynamics  may  play  in  new  fisheries  man¬ 
agement  approaches,  including  the  design  and  implementation  of  marine  protected 
areas  (Palumbi,  2003;  Shanks  et  al.,  2003).  However,  studying  larval  dispersal  (and 
therefore  its  role  in  connectivity)  is  challenging  due  to  small  larval  sizes,  high  dilu¬ 
tion  rates,  and  high  larval  mortality  rates  (Thorson,  1950,  1960).  In  recent  years, 
the  use  of  artificial  and  natural  tags  to  track  marine  larvae  has  been  explored  (e.g. 
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Levin,  1990;  Thorrold  et  al. ,  2002).  One  type  of  natural  tag  that  may  be  useful  for 
identifying  natal  origins  is  the  elemental  signature  recorded  in  biogenic  carbonate. 
This  technique  relies  on  the  observation  that  some  elements  are  incorporated  into 
the  calcium  carbonate  matrix  in  amounts  related  to  the  dissolved  concentrations  or 
physical  properties  of  the  ambient  water  (e.g.  Bath  et  al.,  2000;  Vander  Putten  et  al., 
2000;  Elsdon  and  Gillanders,  2003).  Provided  water  chemistry  or  temperature  is  sig¬ 
nificantly  different  between  natal  habitats,  such  variation  can  serve  as  a  natural  tag, 
or  signature,  of  the  geographic  origin  of  organisms.  The  use  of  geochemical  signatures 
in  fish  otoliths  as  natural  tags  for  population  studies  is  well  established  (Campana 
and  Thorrold,  2001).  Recent  efforts  have  expanded  the  use  of  elemental  tags  to  inver¬ 
tebrates  including  decapods  (DiBacco  and  Levin,  2000),  gastropods  (Zacherl  et  al., 
2003a),  bivalves  (Becker  et  al.,  2005,  2007),  and  cephalopods  (Arkhipkin  et  al.,  2004; 
Zumholz  et  al.,  2007). 

Although  it  is  well  established  that  biogenic  carbonate  composition  reflects  the 
elemental  concentrations  of  waters  where  it  was  formed,  additional  factors  are  known 
to  influence  elemental  uptake  (Dodd,  1967;  Wilbur,  1972).  Environmental  conditions, 
such  as  temperature,  salinity,  oxygen  levels,  and  nutrient  load,  may  alter  the  direct 
relationship  between  environmental  and  carbonate  composition.  Also,  physiological 
and  biochemical  factors  may  affect  the  signature,  including  differences  in  uptake  be¬ 
tween  species  or  between  stages  within  a  species.  These  factors  may  potentially  act  in 
combination,  thereby  complicating  a  simple  relationship  between  ambient  elemental 
concentration  and  carbonate  composition.  Understanding  the  role  each  factor  plays 
in  determining  carbonate  chemistry  is  critical  to  the  appropriate  design  and  inter¬ 
pretation  of  studies  that  use  natural  elemental  signatures  as  a  proxy  for  geographic 
origins  or  to  reconstruct  environmental  conditions. 

Our  ultimate  goal  is  to  develop  elemental  signatures  in  shell  as  tags  of  natal  habitat 
for  the  softshell  clam,  Mya  arenaria ,  a  commercially  important  bivalve  commonly 
found  in  New  England  estuaries.  In  this  study,  we  measured  larval  and  juvenile 
shell  elemental  composition  of  M.  arenaria  reared  in  the  laboratory  under  controlled 
temperature  and  salinity  to  understand  their  effects  on  elemental  uptake  and  therefore 
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allow  us  to  interpret  variability  in  elemental  signatures.  Temperature  and  salinity  are 
likely  to  vary  on  small  time  scales  in  estuaries  (Bowden,  1980),  and  are  known  to 
influence  elemental  uptake  into  biogenic  carbonate. 


First,  we  determined  how  the  relationship  between  elemental  composition  of  sea¬ 
water  and  larval  shell  varied  with  temperature  and  salinity.  Many  studies  have  de¬ 
tailed  the  various  factors  controlling  elemental  compositions  of  fish  otoliths  (e.g.  Hoff 
and  Fuiman,  1995;  Bath  et  al.,  2000;  Arai  et  ah,  2003;  Martin  et  ah,  2004)  and 
adult  mollusc  shells  (e.g.  Rucker  and  Valentine,  1961;  Dodd,  1965;  Lerman,  1965; 
Dodd  and  Crisp,  1982;  Pitts  and  Wallace,  1994;  Lazaret.h  et  ah,  2003;  Zacherl  et  ah, 
2003b).  For  mollusc  larvae,  however,  shell  uptake  has  been  examined  only  in  one 
gastropod  (Zacherl  et  ah,  2003b)  and  in  no  bivalves.  We  explored  the  relationships 
between  temperature  and  salinity  and  uptake  of  Mg,  Mn,  Sr,  Ba,  and  Pb  based  on 
the  demonstrated  utility  of  these  elements  in  previous  uptake  experiments. 


Second,  we  determined  how  the  relationship  between  elemental  composition  of 
seawater  and  juvenile  shell  varied  with  temperature  and  salinity,  and  compared  the 
uptake  into  shells  of  larval  and  juvenile  clams  to  look  for  predictable  correlations. 
The  physiological  aspects  of  shell  formation  in  larval  and  juvenile  bivalves  are  likely 
to  vary  due  to  morphological  differences  in  these  two  life  stages.  Although  the  mantle 
(or  the  precursor  to  the  mantle,  the  shell  gland)  is  responsible  for  shell  production  in 
both  larval  and  post-larval  clams,  there  are  drastic  changes  in  the  tissue  during  devel¬ 
opment  that  might  result  in  differences  in  shell  composition.  The  differences  between 
larval  and  juvenile  shell  composition  may  offer  insight  into  the  physiological  factors 
potentially  affecting  elemental  incorporation  into  shell.  If  temperature  or  salinity 
affects  incorporation  of  elements  into  carbonate,  investigators  may  need  to  consider 
these  factors  when  interpreting  environmental  variables  based  on  shell  chemistry. 
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4.2  Methods 


4.2.1  Clam  rearing 

Adult  Mya  arenaria  were  obtained  from  Cotuit,  Massachusetts  in  April  2006  and 
were  transported  to  the  Environmental  Systems  Laboratory  (ESL)  at  Woods  Hole 
Oceanographic  Institution,  where  they  were  placed  in  mesh  bags  and  suspended  in 
a  750  L  tank  with  filtered  seawater.  Spawning  activity  commenced  approximately 
one  hour  later,  and  the  tank  was  left  undisturbed  to  allow  spawning  to  finish  and 
for  fertilization  to  take  place.  After  five  hours  the  adult  clams  were  removed,  and 
the  tank  contents  were  filtered  through  a  35  //m  nylon  mesh  sieve  to  concentrate  the 
larvae  into  a  small  volume  (20  L).  We  subsampled  and  counted  larvae  to  obtain  a 
total  estimate  of  16  million  trochophore  larvae.  Although  polyspermy  was  evident 
due  to  the  high  ratio  of  sperm  to  eggs,  most  larvae  were  well  formed  and  regular  in 
appearance  and  behavior. 

Trochophore  larvae  were  placed  into  one  of  24  12  L  high-density  polyethylene 
(HDPE)  tanks  at  a  density  of  45  larvae  ml-1.  The  tanks  were  divided  among  six 
treatments  representing  all  combinations  of  two  salinities  (22.5,  30  PSU;  hereafter 
low  and  high)  and  three  temperatures  (15,  20,  and  24  °C),  with  four  replicate  tanks 
per  treatment  combination.  These  temperatures  and  salinities  represent  the  range 
of  conditions  experienced  by  this  species  during  the  summer  months  in  its  natural 
estuarine  habitats.  Treatment  temperatures  were  maintained  by  placing  tanks  in  large 
water  baths.  Water  baths  for  the  15  and  20° C  treatments  received  a  continuous  supply 
of  15  and  20°C  water,  respectively,  regulated  for  the  entire  seawater  supply  line  by 
chillers  and  heaters.  The  water  bath  for  24°C  treatment  tanks  was  maintained  using 
a  120-volt  tank  heater  regulated  by  a  thermostat  (Process  Technology  Easy  Plug™ 
Heater  with  Digital  Controller,  1800W).  Seawater  for  all  tanks  was  obtained  from 
the  in-house  seawater  supply  line,  which  pumps  water  from  Vineyard  Sound  100 
m  offshore  of  ESL  at  a  water  depth  of  4  m.  All  seawater  was  filtered  to  remove 
particles  >  1  /tin  before  use.  Salinity  was  reduced  for  the  low-salinity  treatments 
by  adding  ultrapure  H20  (Milli-Q)  to  Vineyard  Sound  seawater.  Tank  temperatures 
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and  salinities  were  measured  every  other  day,  and  water  samples  were  taken  weekly 
to  determine  ambient  elemental  ratios  over  the  course  of  the  experiment  . 

Subsamples  of  larvae  were  removed  from  each  tank  after  36  h  for  larval  shell 
analysis.  The  remaining  individuals  were  raised  for  60  days  under  treatment  condi¬ 
tions,  with  complete  water  changes  every  two  days.  Temperature  and  conductivity 
measurements  were  taken  daily  for  each  tank.  New  tank  water  was  adjusted  to  the 
appropriate  temperature  and  salinity  before  larvae  were  added.  For  the  first  10  days 
the  larvae  were  allowed  to  swim  freely  in  tanks.  Water  changes  consisted  of  filtering 
out  larvae  through  a  series  of  mesh  sieves,  then  returning  them  to  the  cleaned  tanks. 
Most  larvae  were  large,  negatively  buoyant  pediveligers  in  the  20°C  and  24°C  tanks 
at  8  d  and  in  the  15°C  tanks  at  10  d.  After  10  days,  larvae  were  retained  within 
mesh  sieves  suspended  in  the  tank,  which  restricted  settlement  to  sieve  surfaces  and 
facilitated  tank  cleaning.  Larvae  metamorphosed  into  juveniles  within  two  weeks  of 
spawning.  Larvae  and  juveniles  were  fed  a  mix  of  live  Isochrysis  sp.  and  concent  rated 
algae  (Instant  Algae©  Shellfish  Diet  1800)  one  to  two  times  daily  based  on  clearance 
rates.  All  treatments  and  tanks  received  the  same  food  mixture  over  the  course  of 
the  experiment.  After  60  days,  clams  ranged  in  size  from  0.5  mm  to  10  mm;  they 
were  removed  from  their  tanks  and  frozen  until  the  shells  were  cleaned  and  prepared 
for  analysis. 

4.2.2  Sample  preparation 

Seawater  preparation 

We  monitored  the  elemental  composition  of  ambient  seawater  in  the  tanks  during  the 
experiment.  Water  samples  from  each  tank  were  collected  weekly  over  the  course 
of  the  experiment.  The  samples  were  vacuum  filtered  using  acid-washed  plastic 
Millipore©  funnels  with  0.2  //in  cellulose  nitrate  membrane  filters.  Samples  were 
then  transferred  to  acid- washed  HDPE  bottles,  acidified  to  pH  2  using  ultrapure 
HN03  (Seastar  Chemicals  Inc.),  and  refrigerated  for  up  to  two  months.  Seawater 
samples  were  prepared  for  elemental  analysis  by  diluting  50-fold  with  ultrapure  2% 
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hno3. 


Larval  shell  preparation 

Subsamples  of  roughly  1000  individual  larvae  were  removed  from  each  tank  36  hours 
after  spawning.  This  allowed  for  the  larval  prodissoconch  I  (PI)  shell  to  set  under 
treatment  conditions,  with  minimal  addition  of  post-PI  shell  material.  Individual 
larvae  were  pooled  into  one  vial  per  tank.  Shells  were  soaked  in  1%  NaCIO  (bleach) 
for  6  hours  to  remove  organic  material.  The  cleaning  solution  was  removed  and 
shells  were  transferred  to  clean  0.5  ml  centrifuge  tubes.  Tubes  were  centrifuged  to 
concentrate  larval  shells,  and  then  were  rinsed  three  times  with  ultrapure  H20  (Milli- 
Q).  Shells  were  soaked  in  200  /d  of  1%  ultrapure  H2O2  (ULTREX*  II  Ultrapure 
through  VWR  International)  buffered  to  IN  in  ultrapure  NaOH  (EMD  Chemicals 
through  VWR  International)  for  five  minutes,  then  rinsed  three  times  with  ultrapure 
H2O,  transferred  to  clean  acid-washed  1  ml  vials,  and  then  rinsed  four  more  times. 
Vials  were  transferred  to  a  laminar  flow  hood  to  dry  overnight.  After  all  fluid  had 
evaporated,  1  ml  of  2%  HNO3  was  added  to  each  vial  to  dissolve  the  shells  for  solution- 
based  ICP-MS  analyses. 

Juvenile  shell  preparation 

Juvenile  clams  were  placed  in  individual  acid- washed  vials  using  acid- washed  plastic 
forceps,  and  sonicated  briefly  to  remove  tissue  and  debris.  Individuals  were  rinsed 
three  times  with  ultrapure  H20,  and  then  soaked  for  10  minutes  at  80°C  in  1%  H2O2 
solution  buffered  in  IN  ultrapure  NaOH  to  remove  organic  material.  Afterward,  shells 
were  rinsed  three  times  with  ultrapure  H20,  transferred  to  clean,  acid-washed  vials, 
rinsed  four  times  with  ultrapure  H20,  then  left  to  dry  overnight  under  a  laminar  flow 
hood.  After  the  shells  were  dry,  they  were  dissolved  for  solution-based  analyses.  The 
retained  larval  shell  was  not  removed  prior  to  dissolution  of  the  juvenile  shell  since 
any  contribution  to  elemental  signal  by  the  larval  shell  (<  0.05  //g)  was  assumed  to 
be  negligible  compared  to  the  more  massive  juvenile  material  (mean  weight  =  2.5 
mg).  We  measured  elemental  ratios  of  seawater  from  experimental  tanks  to  calcu- 
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late  discrimination  coefficients  for  both  larval  and  juvenile  Mya  arenaria  (Morse  and 
Bender,  1990).  The  discrimination  coefficient  D  describes  the  relationship  between 
the  elemental  composition  of  carbonate  and  of  its  water  of  formation: 

d  _  (Element:Ca) carbonate 
( Element :  Ca)  seawater 

where  D  <  1  indicates  discrimination  against  inclusion  of  the  element,  D  =  1  suggests 
there  is  no  discrimination,  and  D  >  1  suggests  elemental  enrichment. 

ICP-MS  analyses 

Shell  material  and  water  samples  were  analyzed  using  a  Thermo- Finnigan  MAT  Ele¬ 
ment^  magnetic  sector  field  inductively  coupled  plasma  mass  spectrometer  (ICP-MS). 
To  correct  for  mass  bias  and  instrument  drift,  standards  were  run  periodically  during 
the  analyses.  For  shell  analyses,  the  solution-based  standards  included  an  aragonitic 
otolith  reference  material  (Yoshinaga  et  ah,  2000)  and  the  certified  reference  material 
(CRM)  FEBS-1  (Sturgeon  et  ah,  2005).  For  seawater  analyses,  CASS-4  Nearshore 
Seawater  Reference  Material  was  used  as  the  standard  (National  Research  Council 
Canada  CRM). 

During  analyses  we  monitored  25Mg,  55Mn.  88Sr,  114Cd,  138Ba,  and  208Pb  in  low 
resolution  mode.  These  elements  were  chosen  because  of  their  previous  use  in  suc¬ 
cessful  elemental  tagging  studies  (Kalish,  1989:  Elsdon  and  Gillanders,  2002:  Zacherl 
et.  ah,  2003b;  Martin  and  Thorrold,  2005).  Molar  ratios  of  each  element  to  Ca  (here¬ 
after  E:Ca)  were  calculated  using  mass  bias  corrections  calculated  from  calibration 
standards.  Limits  of  detection  (LOD)  were  calculated  as  the  ratio  of  three  standard 
deviations  of  the  blank  intensity  to  the  average  blank-subtracted  sample  intensity. 
For  seawater  and  solution-based  analyses  respectively,  sample  intensities  were  as  fol¬ 
lows:  >  10000  and  >  1200  times  the  LOD  for  25Mg;  >  700  and  >  10000X  LOD  for 
48Ca;  >16  and  >  57 x  LOD  for  55Mn;  >  5500  and  >  llOOx  LOD  for  88Sr;  >  1600 
and  >  130x  LOD  for  138Ba;  and  >  75  and  >  15x  LOD  for  208Pb.  The  intensities  of 
U4Cd  were  not  routinely  above  detection  limits  and  this  element  was  excluded  from 
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the  remaining  analyses.  Elemental  ratios  that  were  more  than  two  standard  devia¬ 
tions  from  the  global  mean  were  eliminated  from  further  analysis  as  the  anomalous 
values  likely  arose  from  sample  contamination.  This  resulted  in  removal  of  6.5%  (66 
of  1010)  of  all  shell  measurements  and  7.4%  (74  of  1001)  of  all  seawater  measure¬ 
ments.  Measured  precision  (%  relative  standard  deviation)  of  E:Ca  during  analysis  of 
juvenile  shells  (n  =  32),  larval  shells  (n  =  10),  and  seawater  (n  =  52)  respectively  was 
2.4%,  4.8%,  and  7.9%  for  Mg:Ca;  8.6%,  11.9%,  and  7.1%  for  Mn:Ca;  0.31%,  0.12%, 
and  0.47%  for  Sr:Ca;  1.8%,  1.4%,  and  2.1%  for  Ba:Ca;  and  2.5%,  7.1%,  and  3.5%  for 
Pb:Ca. 

4.2.3  Statistical  analyses 

Although  seawater  for  all  tanks  was  obtained  from  the  same  supply  line,  it  was  impor¬ 
tant  to  quantify  any  differences  among  tanks  as  these  differences  may  have  influenced 
patterns  in  shell  uptake.  We  used  two-way  analysis  of  variance  (ANOVA)  to  test 
for  the  effects  of  temperature  and  salinity  on  seawater  E:Ca,  with  individual  tank 
analyses  for  each  week  as  the  unit  of  measurement  (n  =  216).  Both  temperature  and 
salinity  were  treated  as  fixed  variables. 

For  each  shell  measurement,  E:Ca  was  translated  into  a  discrimination  coefficient. 
Larval  shell  discrimination  coefficients  were  calculated  by  dividing  E:Ca  in  the  shell  by 
E:Ca  for  corresponding  tank  seawater  sampled  on  day  2  of  the  experiment.  Juvenile 
shell  discrimination  coefficients  were  calculated  by  dividing  the  shell  E:Ca  by  the 
corresponding  mean  tank  seawater  E:Ca  (averaged  from  samples  taken  over  the  nine 
weeks  of  the  experiment).  We  used  discrimination  coefficients  rather  than  elemental 
ratios  of  shell  for  statistical  analyses  because  the  discrimination  coefficient  accounts 
for  possible  differences  in  water  chemistry  among  tanks.  We  used  a  three-way  ANOVA 
to  explore  the  effects  of  temperature,  salinity  and  stage  (i.e.  larval  versus  juvenile 
shell)  on  uptake.  There  were  significant  stage  effects  on  discrimination  coefficients 
(p  <  0.01  for  Mg;  p  <  0.0001  for  Mn,  Sr,  Ba,  Pb;  data  not  shown)  and  consequently 
we  chose  to  analyze  larval  and  juvenile  shells  separately  using  two-way  A NOVAs  to 
test  for  temperature  and  salinity  effects.  Finally,  we  tested  for  differences  between 
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Figure  4.1:  Seawater  elemental  ratios  (mean  of  replicate  tanks  ±SE)  for  treatments 
over  experimental  week  for  (a)Mg:Ca,  (b)Mn:Ca,  (c)Sr:Ca,  (d)Ba:Ca,  and  (e)Pb:Ca. 
Gray,  dotted,  and  black  lines  are  15°C,  20°C,  and  24°C  treatments  respectively.  Open 
symbols  are  low-salinity  treatments,  and  closed  symbols  are  high-salinity  treatments. 


larval  and  juvenile  uptake  using  simple  linear  correlation  analysis  of  discrimination 
coefficients. 


4.3  Results 

4.3.1  Rearing  conditions 

Although  wre  made  efforts  to  optimize  conditions  for  survival  and  growth,  survivor¬ 
ship  was  highly  variable  between  tanks.  Some  tanks  had  hundreds  of  surviving  clams 
after  the  60  day  interval  of  the  experiment  ,  while  two  tanks  had  no  survivors  (Table 
4.1).  We  therefore  analyzed  juvenile  clam  shells  from  22  of  24  tanks,  with  at  least 
three  replicate  tanks  for  each  treatment.  Elemental  ratios  were  obtained  from  a  total 
of  23  larval  shell  samples  (from  23  tanks)  and  174  juvenile  shells  (from  22  tanks)  to 
examine  temperature,  salinity,  and  stage  effects  on  uptake.  Seawater  elemental  com¬ 
position  varied  among  treatments  and  tanks  (Fig.  4.1).  A  two-way  ANOVA  showed 
that  elemental  ratios  significantly  varied  with  salinity  for  Mg  and  Mn  (low  salinity  > 
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Tank 

N 

Temp 

Sal 

Mg:Ca 

Mn:Ca 

Sr:Ca 

Ba:Ca 

Pb:Ca 

1 

32 

15.4  (0.10) 

22.6  (0.06) 

5.00  (0.03) 

30.3  (4.74) 

8.53  (0.02) 

14.7  (0.83) 

9.17(1.03) 

2 

2 

15.5  (0.10) 

22.6  (0.07) 

5.03  (0.02) 

27.5  (3.80) 

8.52  (0.02) 

16.9(1.54) 

15.9  (3.40) 

3 

30 

15.3  (0.09) 

22.5  (0.06) 

5.05  (0.02) 

21.4  (5.01) 

8.54  (0.02) 

14.4(1.09) 

8.61  (1.46) 

4 

0 

15.3  (0.10) 

22.5  (0.05) 

5.02  (0.02) 

18.7(6.34) 

8.54  (0.03) 

16.6(1.58) 

14.1  (5.02) 

5 

5 

15.3  (0.09) 

30.1  (0.04) 

5.00  (0.02) 

20.0  (5.20) 

8.53  (0.02) 

15.0  (0.77) 

10.0  (2.37) 

6 

5 

15.4  (0.09) 

30.0  (0.04) 

5.03  (0.02) 

15.6  (3.76) 

8.56  (0.02) 

13.2  (0.85) 

6.29(1.10) 

7 

6 

15.4  (0.09) 

30.0  (0.05) 

5.03  (0.02) 

17.9  (4.09) 

8.55  (0.03) 

13.0(0.61) 

11.6  (3.74) 

8 

7 

15.4  (0.09) 

30.0  (0.05) 

5.05  (0.01) 

18.5(3.58) 

8.51  (0.02) 

15.2  (0.79) 

12.2  (3.19) 

9 

1 

19.4  (0.04) 

22.3  (0.07) 

5.04  (0.01) 

17.9  (3.27) 

8.56  (0.02) 

15.1  (1.06) 

10.0(1.12) 

10 

>100 

19.4  (0.04) 

22.4  (0.07) 

5.09  (0.04) 

18.6(4.94) 

8.56(0.03) 

14.9(1.38) 

15.3  (4.12) 

11 

>100 

19.4  (0.04) 

22.4  (0.07) 

5.09  (0.03) 

13.6(3.55) 

8.64  (0.03) 

15.1  (1.41) 

9.37  (2.27) 

12 

>100 

19.5  (0.04) 

22.3  (0.05) 

5.05  (0.02) 

14.7  (3.98) 

8.57  (0.03) 

15.1  (1.19) 

9.32  (1.53) 

13 

>100 

19.4  (0.03) 

29.7  (0.04) 

5.02  (0.02) 

13.6  (3.00) 

8.52  (0.01) 

14.9  (0.44) 

6.80(1.75) 

14 

4 

19.4  (0.04) 

29.8  (0.05) 

5.00  (0.02) 

15.8  (4.09) 

8.51  (0.01) 

15.2(1.08) 

9.12(2.01) 

15 

8 

19.4  (0.03) 

29.8  (0.06) 

5.01  (0.02) 

19.6  (3.45) 

8.51  (0.01) 

18.3  (2.58) 

11.2  (2.09) 

16 

4 

19.4  (0.03) 

29.8  (0.06) 

5.02  (0.01) 

14.9  (4.13) 

8.53  (0.01) 

15.2  (0.91) 

12.0(4.03) 

17 

0 

23.9  (0.05) 

22.5  (0.10) 

5.06  (0.01) 

13.8(3.13) 

8.55  (0.01) 

17.5  (2.21) 

12.9  (2.29) 

18 

15 

23.9  (0.05) 

22.7  (0.27) 

5.06  (0.01) 

16.2  (4.11) 

8.60  (0.01) 

16.3  (1.32) 

9.99  (1.52) 

19 

1 

23.8  (0.05) 

22.4  (0.09) 

5.05  (0.01) 

15.4  (3.13) 

8.55  (0.01) 

15.0  (0.77) 

14.1  (3.46) 

20 

>100 

23.9  (0.05) 

22.4  (0.09) 

5.05  (0.02) 

15.5  (4.25) 

8.61  (0.01) 

16.8(1.75) 

11.5  (2.26) 

21 

>100 

23.9  (0.05) 

29.5  (0.25) 

5.03  (0.01) 

12.8(3.26) 

8.58  (0.03) 

14.9  (0.68) 

1 1.0  (2.04) 

22 

>100 

23.8  (0.05) 

29.7  (0.06) 

5.03  (0.02) 

11.6  (3.66) 

8.61  (0.03) 

14.5  (0.17) 

8.31  (1.91) 

23 

40 

24.0  (0.05) 

29.7  (0.08) 

5.01  (0.03) 

12.6(2.65) 

8.61  (0.04) 

15.6(1.03) 

9.45  (2.26) 

24 

>100 

23.9  (0.05) 

29.5  (0.25) 

5.02  (0.01) 

12.2  (3.30) 

8.56  (0.01) 

15.5(0.38) 

18.7(5.83) 

Table  4.1:  Number  of  surviving  clams  (N)  at  the  end  of  the  60  day  experiment  and 
mean  water  temperature  (Temp.,  °C  [±SE],  n  =  30),  salinity  (Sal.,  PSU  [±SE], 
n  =  30),  dissolved  Mg:Ca  (mol  mol-1,  ±SE),  Mn:Ca  (//mol  mol-1,  ±SE),  Sr:Ca 
(mmol  mol-1,  ±SE),  Ba:Ca  (/rmol  mol-1,  ±SE),  and  Pb:Ca  (nmol  mol-1,  ±SE)  for 
treatment  tanks  (n  =  9  for  each  tank). 
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Source 

DF 

Mg:Ca 

MS  F 

Mn:Ca 

MS  F 

Sr:Ca 

MS  F 

Ba:Ca 

MS  F 

Pb.Ca 

MS  F 

T 

2 

7  xlO'3 

2.04 

968 

7.62** 

4  x  1 0 2 

10.7** 

11.1 

0.903 

29.4 

0.46 

S 

1 

4  xlO 2 

12.3** 

510 

4.02* 

1  xlO2 

2.98 

28.5 

2.33 

92.1 

1.44 

Tx  S 

2 

8  x  1 0 1 

2.42 

164 

1.29 

3  xlO'2 

7.02** 

21.0 

1.72 

5.51 

0.09 

Error 

195 

3  xlO 3 

127 

4  x  1 0 3 

12.3 

63.9 

Table  4.2:  Results  of  ANOVA  testing  the  effects  of  temperature  (T)  and  salinity  (S) 
on  tank  water  sample  dissolved  elemental  ratios.  All  tank  measurements  for  all  weeks 
were  used  as  the  experimental  unit  (n  =  216).  F  statistics  are  significant  at  the  level 
of  *p  <  0.01  or  **p  <  0.001. 

high  salinity),  and  with  temperature  for  Mn  and  Sr  (Table  4.2).  Seawater  elemental 
composition  also  varied  temporally  (Fig.  4.1).  In  order  to  account  for  this  variability 
in  seawater  chemistry  across  treatments,  we  calculated  discrimination  coefficients  for 
shells  using  seawater  values  from  week  one  only  for  larval  discrimination  coefficients 
(Table  4.3),  and  using  seawater  values  averaged  over  the  nine-week  experiment  for 
juvenile  discrimination  coefficients  (Table  4.1).  Whole-shell  solution  analysis  of  juve¬ 
niles  resulted  in  shell  E:Ca  values  that  represent  an  average  of  elemental  uptake  over 
the  experiment.  By  using  the  nine-week  average  for  each  tank  to  calculate  juvenile 
shell  discrimination  coefficients,  we  are  accounting  for  this  averaging. 


4.3.2  Temperature  &;  salinity  effects  on  uptake 

Larval  shell 

Across  all  treatments  and  tanks,  the  average  (±SE)  larval  shell  elemental  ratios  were 
3.1  ±0.83  mmol  mol-1  for  Mg:Ca,  44.8  ±4.1  //mol  mol-1  for  Mn:Ca,  2.5  ±0.04  mmol 
mol-1  for  Sr:Ca,  3.37  ±  0.30  //mol  mol-1  for  Ba:Ca,  and  157  ±  17.6  nmol  mol-1  for 
Pb:Ca.  Average  discrimination  coefficients  were  less  than  one  for  Mg  (6.2  x  10-4  ± 
1.7  x  10~4),  Sr  (0.30  ±  0.005),  and  Ba  (0.24  ±  0.02),  and  greater  than  one  for  Mn 
(1.86  ±  0.19)  and  Pb  (22.3  ±  3.8). 

Temperature  significantly  affected  larval  uptake  of  Ba  and  Mn  (Table  4.4;  Fig.  4.2). 
Mean  Defl  for  larval  shells  from  15°C  tanks  was  higher  than  those  from  20°C  or  24°C 
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Tank 

Temp 

Sal 

Mg:Ca 

Mn:Ca 

Sr:Ca 

Ba:Ca 

Pb:Ca 

i 

14.9  (0.15) 

22.7  (0.06) 

4.91 

44.0 

8.53 

17.49 

11.4 

2 

15.1  (0.20) 

22.8  (0.45) 

5.10 

37.77 

8.45 

15.42 

17.1 

3 

14.9  (0.12) 

22.9  (0.34) 

5.15 

31.40 

8.46 

12.18 

9.88 

4 

14.8(0.12) 

22.6  (0.30) 

5.06 

32.47 

8.60 

12.19 

10.5 

5 

15.0  (0.23) 

30.1  (0.11) 

5.10 

17.68 

8.46 

13.24 

5.20 

6 

15.2  (0.12) 

30.0  (0.10) 

5.05 

17.15 

8.55 

13.57 

5.74 

7 

15.1  (0.07) 

30.2  (0.26) 

4.86 

16.41 

8.44 

12.98 

6.83 

8 

15.0  (0.09) 

30.0  (0.22) 

5.11 

19.78 

8.56 

14.72 

7.43 

9 

19.5  (0.26) 

22.1  (0.11) 

5.09 

23.65 

8.53 

12.59 

6.27 

10 

19.5  (0.17) 

22.1  (0.15) 

5.06 

24.85 

8.60 

20.97 

4.92 

11 

19.5  (0.26) 

22.2  (0.06) 

5.04 

24.37 

8.51 

14.54 

7.31 

12 

19.5  (0.26) 

22.0  (0.08) 

5.09 

34.48 

8.63 

1205 

3.93 

13 

19.5  (0.20) 

29.7  (0.17) 

5.05 

21.74 

8.60 

14.62 

7.95 

14 

19.5  (0.21) 

29.4  (0.05) 

5.01 

17.09 

8.55 

13.71 

6.25 

15 

19.5  (0.19) 

29.3  (0.09) 

5.06 

15.92 

8.58 

13.97 

3.98 

16 

19.5  (0.19) 

29.4  (0.19) 

4.99 

24.56 

8.54 

15.37 

33.5 

17 

23.8  (0.38) 

22.3  (0.30) 

5.09 

22.91 

8.57 

18.91 

4.06 

18 

23.8(0.19) 

22.1  (0.26) 

4.99 

31.95 

8.48 

15.61 

27.5 

19 

23.4  (0.20) 

22.2  (0.44) 

5.01 

19.09 

8.69 

13.08 

8.11 

20 

23.5  (0.13) 

22.1  (0.19) 

5.03 

27.28 

8.53 

18.17 

15.0 

21 

23.6(0.15) 

29.4  (0.08) 

4.99 

23.08 

8.50 

13.53 

5.12 

22 

23.5  (0.09) 

29.7  (0.31) 

5.07 

27.23 

8.52 

13.51 

3.56 

23 

23.7  (0.15) 

29.1  (0.34) 

4.98 

25.51 

8.53 

13.45 

6.17 

24 

23.8  (0.35) 

29.7  (0.09) 

4.96 

24.84 

8.46 

14.87 

7.33 

Table  4.3:  Mean  water  temperature  (Temp.,  °C  [±SE],  n  =  3),  mean  salinity  (Sal., 
PSU  [±SE],  n  =  3),  and  dissolved  Mg:Ca  (mol  mol-1),  Mn:Ca  (//mol  mol-1),  Sr:Ca 
(mmol  mol-1),  Ba:Ca  (//mol  mol-1),  and  Pb:Ca  (nmol  mol-1)  for  treatment  tanks 
during  the  first  week  only. 
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Source 

DF 

MS 

F 

MS 

F 

MS 

F 

MS 

F 

MS 

F 

Larval  shell 

T 

2 

7.3xl0'7 

1.34 

2.50 

9.50* 

9.9x1 0-4 

3.45 

0.081 

27.7** 

357 

1.03 

s 

1 

4.4x1  O'7 

0.81 

2.59 

9.85* 

1.0x1 0'3 

3.75 

9.6x10-' 

3.26 

537 

1.55 

Tx  S 

2 

2.5x1  O'7 

0.46 

2.19 

8.33* 

1.2x10 1 

4.06 

7.7x10-' 

2.62 

184 

0.53 

Error 

15 

5.5x1 0-7 

0.26 

2.9x1 0-4 

2.9x10-' 

347 

Juvenile  shell 

T 

2 

2.3x10-'° 

0.82 

0.39 

6.43* 

3.4x1  O'4 

4.55 

6.6x10-' 

9.24* 

0.72 

10.5* 

S 

1 

6.4x10" 

0.23 

1 .0x10 3 

0.01 

3.9x10-' 

52.2** 

1.0x1  O'6 

0.001 

0.01 

0.07 

Tx  S 

2 

0.54 

0.54 

0.19 

3.14 

6.5x1  O'4 

8.64* 

8.3x1  O'4 

1.17 

0.04 

0.58 

Error 

15 

2.8x1 0'° 

0.06 

7.5x10' 

7.1X10-4 

0.68 

Table  4.4:  Results  of  ANOVA  testing  the  effects  of  temperature  (T)  and  salinity  (S) 
on  larval  and  juvenile  discrimination  coefficients.  F  statistics  are  significant  at  the 
level  of  *p  <  0.01  or  **p  <  0.001. 


Figure  4.2:  Mean  discrimination  coefficients  (±SE)  for  larval  shell  samples,  averaged 
over  replicate  tanks  within  treatments  for  the  five  elements,  (a)  Da./9,  (b)  D A/n  (c) 
Dsr,  (d)  D«a,  and  (e)  D/^.  Open  circles  are  low-salinity  treatments  and  closed  circles 
are  high-salinity  treatments. 
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tanks  while  shells  from  15°C  tanks  had  higher  mean  Da/„  than  both  20°  C  and  24°  C 
tanks.  Salinity  significantly  affected  Mn  uptake  in  larval  shell:  high-salinity  tanks  had 
significantly  higher  D Mn  than  low-salinity  tanks.  There  was  a  significant  interaction 
between  salinity  and  temperature  in  larval  shells,  caused  by  a  variable  effect  of  salinity 
depending  on  temperature:  low-salinity  tanks  had  lower  T>Mn  than  high-salinity  tanks 
at  15°C  and  24°C,  while  this  trend  was  reversed  at  20°C. 

Juvenile  shell 

The  average  (±SE)  juvenile  shell  elemental  ratios  across  all  treatments  and  tanks  were 
0.67  ±  0.02  mmol  mol-1  for  Mg:Ca,  14.3  ±  0.92  pmol  mol-1  for  Mn:Ca,  1.99  ±  0.04 
mmol  mol-1  for  Sr:Ca,  1.52  ±0.11  /rmol  mol-1  for  Ba:Ca,  and  4.56  ±0.65  nmol  mol-1 
for  Pb:Ca.  All  of  the  average  discrimination  coefficients  for  juvenile  shell  were  less 
than  one  (1.3  x  10~4±1.0  x  10-5  for  Da 0.88  ±0.13  for  D 0.23  ±0.003  for  Dsr, 
0.10  ±  0.008  for  DSo,  and  0.45  ±  0.13  for  Dp*,). 

Temperature  significantly  affected  uptake  of  Mn,  Ba  and  Pb  in  juvenile  shell 
(Table  4.4:  Fig.  4.3).  Tukey-Kramer  pairwise  comparisons  indicated  the  significant 
relationships  were  as  follows:  mean  Da/«  was  higher  for  individuals  from  20°C  tanks 
than  from  either  15°C  or  24°C  tanks;  mean  Dpa  was  higher  in  15°C  tanks  than  in  24°C 
tanks;  and  mean  DPl)  was  higher  for  individuals  from  15°C  tanks  than  from  either 
20°C  or  24°C  tanks.  Mean  Dsr  did  not  vary  consistently  with  temperature,  but  was 
significantly  greater  in  low  than  high  salinity.  There  was  a  significant  interaction 
between  salinity  and  temperature  for  Dpr,  likely  caused  by  the  less  distinct  salinity 
effect  at  20°C  than  at  15°C  or  24°C. 

Tests  for  correlations  between  larval  and  juvenile  uptake  in  each  tank  showed 
consistently  higher  discrimination  coefficients  in  larval  than  juvenile  shells  (Fig.  4.4). 
There  was,  however,  little  correspondence  between  larval  and  juvenile  discrimination 
coefficients  between  tanks  for  any  element  except  Ba  ( p  <  0.05,  r  =  0.56)  (Table 
4.5).  The  discrimination  coefficient  used  for  each  tank  was  a  single  value  for  larvae 
(measured  from  shells  pooled  into  a  vial),  and  a  single  or  average  value  for  juveniles 
(calculated  from  1  to  15  shells).  Only  tanks  with  values  for  both  larvae  and  juveniles 
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1.6 


Figure  4.3:  Mean  discrimination  coefficients  (±SE)  for  juvenile  shells,  averaged  over 
replicate  tanks  within  treatments  for  the  five  elements,  (a)  D Mg,  (b)  DMn  (c)  Dsr, 
(d)  Dpa,  and  (e)  Dp*,.  Open  circles  are  low-salinity  treatments  and  closed  circles  are 
high-salinity  treatments. 


r 

P 

^Mg 

-0.166 

0.51 

DMn 

-0.064 

0.80 

Dsr 

0.120 

0.63 

DBa 

0.560* 

0.02 

DPb 

0.164 

0.50 

Table  4.5:  Results  of  analyses  testing  for  correlations  between  juvenile  and  larval 
discrimination  coefficients.  Correlation  coefficients  are  significant  at  the  level  of  *p  < 


0.05. 
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0.28 


Dga  Juveniles  Dpb  Juveniles 


Figure  4.4:  Correlations  comparing  larval  and  juvenile  discrimination  coefficients 
for  (a)  Da ,lg,  (b)  DMn  (c)  D5r,  (d)  DBa,  and  (e)  DPb.  Gray  lines  are  1:1  lines  for 

D  larvae  •  D juveniles  • 

were  used  in  the  correlations  (n  =  18, 19, 19, 18, 19  for  DMg,  DMn,  Dsr,  DBa,  and  Dp*, 
respectively). 

4.4  Discussion 

We  examined  the  effects  of  temperature  and  salinity  on  discrimination  coefficients 
(and  therefore  uptake)  in  larval  and  juvenile  bivalve  shell  in  order  to  improve  our 
ability  to  interpret  natural  variability  in  elemental  signatures.  First,  we  explored 
relations  between  temperature,  salinity,  and  elemental  composition  of  larval  Mya 
arenaria  shells.  We  found  that  uptake  of  Ba  was  affected  by  temperature,  and  uptake 
of  Mn  was  affected  by  temperature  and  salinity.  Second,  we  compared  uptake  in 
larval  and  juvenile  shell  at  various  temperatures  and  salinities.  We  found  significant 
differences  between  larval  and  juvenile  discrimination  coefficients  for  all  five  of  the 
elements  studied.  Correlation  analysis  further  indicated  that  there  was  no  clear  and 
predictable  relationship  between  larval  and  juvenile  uptake  for  any  element  except 
Ba.  Although  we  expected  that  discrimination  coefficients  would  differ  between  larval 
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and  juvenile  shells,  we  were  surprised  to  find  that  temperature  and  salinity  influenced 
discrimination  in  shells  of  the  two  stages  differently. 

4.4.1  Temperature  &  Salinity  Effects  on  Uptake 

Magnesium  Uptake 

Our  observation  that  neither  temperature  nor  salinity  had  a  significant  effect  on  Dm9 
in  larval  or  juvenile  My  a  arenaria  shells  is  in  contrast  to  results  from  some  previous 
studies  of  bivalve  calcite,  (Rucker  and  Valentine,  1961;  Dodd,  1965)  in  which  a  positive 
correlation  to  temperature  was  reported.  Lerman  (1965),  however,  reported  no  effect 
of  temperature  on  D Mg  in  oyster  shells.  We  found  no  prior  studies  of  salinity  effects 
on  D Mg  in  bivalve  shells.  Calcite  and  aragonite  formation  tend  to  differ  with  relation 
to  solution  chemistry  due  to  their  different  crystal  structures,  even  within  the  same 
organism.  Lorens  and  Bender  (1980)  hypothesized  that  Mg  was  regulated  by  mantle 
cells  during  calcite  formation  but  not  during  aragonite  precipitation  in  the  shell  of 
the  bivalve  Mytilus  edulis.  The  radius  of  Mg2+  ions  is  smaller  than  Ca2+  ions,  and 
therefore  alt  hough  Mg  inhibits  formation  of  calcite  due  to  substitution  for  Ca2+  in  the 
crystal  lattice  (Berner,  1975),  it  is  unlikely  to  be  substituted  into  aragonite  (Onuma 
et  al.,  1979;  Speer,  1983).  However  the  extremely  low  discrimination  coefficients  we 
found  suggest  that,  the  amount  of  Mg  in  shell  carbonate  is  highly  regulated  by  the 
organism  during  shell  formation.  Although  the  mechanisms  regulating  discrimination 
of  Mg  during  shell  formation  in  M.  arenaria  remain  are  unresolved,  temporal  changes 
in  temperature  and  salinity  are  unlikely  to  complicate  interpretation  of  Mg  signatures. 

Manganese  Uptake 

The  influence  of  salinity  on  Mn  upt  ake  between  larval  and  juvenile  shells  suggests  that 
Mil  regulation  changes  during  ontogeny.  Ontogenetic  changes  in  discrimination  of  Mn 
have  not  been  reported  in  other  studies  of  biogenic  carbonate,  but  investigations  of 
otolith  chemistry  have  found  complex  relations  between  Mn  uptake  and  temperature 
and  salinity.  Elsdon  and  Gillanders  (2003)  found  no  significant  relationship  between 
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Acanthopagrus  butcheri  otolith  Mn:Ca  and  the  Mn  content  of  the  water  of  formation. 
Martin  and  Thorrold  (2005)  reported  complex  temperature  and  salinity  interactions 
for  Mn:Ca  in  Leiostromus  xanthurus  otoliths.  The  authors  attributed  this  result  to 
biological  processes  that  might  be  mediating  the  amount  of  biologically  available  Mn 
in  the  seawater  of  formation.  For  instance,  microbial  activity  can  contribute  to  the 
formation  of  Mn  oxides  (Sunda  and  Huntsman,  1987;  Klinkhammer  and  McManus, 
2001),  which  removes  Mn  ions  from  solution  and  makes  relationships  between  Mn:Ca 
in  shell  and  seawater  unpredictable. 

Despite  the  unpredictable  nature  of  Mn  uptake,  Mn  has  repeatedly  proven  useful 
for  elemental  tagging  studies  (e.g.  Thorrold  et  al.,  1998,  2001;  Patterson  et.  al.,  2005; 
Becker  et  ah,  2007).  Our  study  suggests  that  experimenters  should  proceed  with 
caution  when  including  Mn  in  the  suite  of  elements  for  tagging  studies  due  to  the 
interactions  among  temperature,  salinity,  and  biology  that  contribute  to  Mn  content 
in  aragonite.  Further,  salinity  effects  on  Mn  discrimination  in  larval  shells  of  estuarine 
species  (where  salinity  varies  on  small  spatial  scales)  may  result  in  variable  signatures 
for  larvae  within  a  particular  site  and  consequent  misclassification  of  individuals  in 
natural  tagging  studies. 

Strontium  Uptake 

Our  observation  that  temperature  did  not  significantly  affect  uptake  of  Sr  into  larval  or 
juvenile  My  a  arenaria  shell  contrasts  with  early  molluscan  studies  generally  reporting 
negative  correlations  between  temperature  and  Sr  (Dodd,  1965;  Hallam  and  Price, 
1968;  Thorn  et  ah,  1995).  However  shells  analyzed  in  these  studies  were  collected 
from  the  field,  and  the  authors  did  not  account  for  differences  in  growth  rates  or 
metabolism  in  their  analyses.  As  a  result,  temperature-associated  variation  in  growth 
rates  or  metabolism  could  have  been  responsible  for  the  Sr  uptake  patterns  rather 
than  temperature  itself.  A  more  recent  study  of  temperature  effects  on  Sr  in  statoliths 
of  the  laboratory-reared  cuttlefish  Sepia  officinalis  that  controlled  for  growth  rates 
also  reported  no  correlation,  suggesting  biological  regulation  of  Sr  ions  outweighs  the 
kinetic  effects  of  temperature  on  Sr  incorporation  (Zumholz  et  ah,  2007).  Zacherl 
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et  al.  (2003b)  came  to  similar  conclusions  in  their  study  of  the  gastropod  Kelletia 
kelletii :  they  found  that  Sr:Ca  correlated  positively  with  temperature  in  protoconchs 
and  negatively  in  larval  statoliths.  They  concluded  that  there  are  strong  biological 
controls  on  Sr  incorporation  into  aragonite,  likely  due  to  the  differential  discrimination 
of  Sr  across  different  biological  membranes  responsible  for  statolith  and  protoconch 
production. 

There  is  less  ambiguity  in  the  literature  for  the  relationship  between  Sr  uptake 
into  molluscan  calcium  carbonate  and  salinity.  The  most  commonly  reported  result 
is  a  lack  of  correlation  (Dodd,  1965;  Hallam  and  Price,  1968;  Stecher  et  al.,  1996; 
Zumholz  et  al.,  2007),  as  we  found  in  our  larval  shells.  Because  the  low- salinity 
treatment  was  achieved  by  diluting  seawater  with  ultrapure  water,  we  did  not  expect 
to  detect  differences  in  Sr:Ca  among  treatments.  However  we  found  Dgr  significantly 
differed  with  salinity  for  juvenile  shells,  which  suggests  that  Sr  concentration,  as  well 
as  Sr:Ca  ratios,  in  ambient  water  influenced  uptake  of  Sr  into  juvenile  shells.  A  similar 
result  was  found  in  otoliths  by  Martin  and  Thorrold  (2005). 

Our  data  suggest  that  molluscan  Sr:Ca  ratios  in  aragonite  of  shells  are  controlled 
by  physiology.  However,  the  exact  mechanism  by  which  these  “vital  effects”  influence 
shell  chemistry  remains  unknown.  There  is  evidence  in  the  literature  that  Sr  up¬ 
take  is  dependent  upon  calcification  rates  in  aragonitic  mollusc  shells  (Purton  et  al., 
1999;  Takesue  and  van  Geen,  2004).  Any  factor  that  affects  calcification  rates  could, 
therefore,  influence  the  interpretation  of  trends  in  Sr:Ca  in  shell  (Likins  et  al.,  1963: 
Stecher  et  al.,  1996;  Carre  et  al.,  2006).  This  may  explain  reports  that  temperature, 
which  is  often  positively  correlated  with  growth  rate,  influences  Sr  uptake.  Gillikin 
et  al.  (2005),  however,  carefully  documented  the  relationship  between  environmental 
factors  and  Sr  uptake  in  two  species  of  marine  bivalves  with  aragonitic  shells,  and 
concluded  that  Sr:Ca  ratios  in  shell  are  regulated  by  biological  processes  rather  than 
thermodynamics.  Studies  of  growth  rate  effects  on  Sr  incorporation  into  Mya  are- 
naria  juvenile  shell  have  produced  both  positive  and  negative  correlations  to  growth 
rate  for  juvenile  shells  reared  in  the  same  conditions  (Strasser  et  al.,  2008).  These 
ambiguous  results  lend  support  to  the  idea  that  Sr  incorporation  may  involve  a  num- 
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ber  of  physical  and  biological  processes.  One  example  of  a  process  that  might  be 
influencing  Sr  incorporation  is  entrapment  of  material  on  the  surface  of  the  crystal 
during  carbonate  formation.  Recent  studies  suggests  that  elemental  fractionation  in 
inorganic  calcite  (Watson,  2004)  and  aragonite  (Gaetani  and  Cohen,  2006)  is  domi¬ 
nated  by  non-equilibrium  processes  including  surface  entrapment.  Indeed,  with  more 
study,  surface  entrapment  may  provide  a.  powerful  unifying  mechanism  for  explain¬ 
ing  apparently  divergent  behavior  of  Sr:Ca  in  fish  otoliths,  bivalve  shells  and  coral 
skeletons. 

Barium  and  Lead  Uptake 

We  found  evidence  for  a  negative  correlation  between  temperature  and  Ba  uptake 
in  both  larval  and  juvenile  shells,  which  is  consistent  with  results  of  Zacherl  et.  al. 
(2003b)  and  Zumholz  et  al.  (2007),  who  both  found  a  negative  correlation  in  the 
aragonite  structures  of  laboratory-reared  molluscans.  We  found  no  significant  rela¬ 
tionship  between  D#0  and  salinity,  which  is  also  consistent  with  results  from  Zumholz 
et  al.  (2007).  Ba  is  likely  to  be  a  useful  addition  to  elemental  tagging  studies,  even  in 
situations  where  its  concentration  does  not  vary  spatially,  as  uptake  into  molluscan 
shell  appears  to  have  a  consistent  relationship  with  temperature. 

Larval  shell  Dpj,  did  not  vary  significantly  with  temperature,  although  juvenile 
shell  Dpt  showed  evidence  of  a  significant  negative  correlation.  Previous  research 
generally  has  shown  that  clams  incorporate  Pb  into  their  shells  in  relation  to  ambient 
water  concentrations  (Babukutty  and  Chacko,  1992;  Pitts  and  Wallace,  1994;  Almeida 
et  al.,  1998;  Boisson  et  al.,  1998),  however  these  studies  do  not  examine  the  relation¬ 
ship  at  different  temperatures  or  salinities.  More  recently  Mubiana  and  Blust  (2007) 
showed  that  Pb  uptake  into  soft  tissue  is  positively  correlated  to  temperature,  but 
they  did  not  examine  the  effects  of  temperature  on  Pb  uptake  into  shell.  It  is  difficult 
to  discern  the  underlying  causes  of  the  significant  relationship  between  temperature 
and  juvenile  shell  DPb  that  we  found.  Therefore  studies  that  use  trace  metals  as  part 
of  unique  geochemical  signatures  should  be  cautious  until  there  is  more  information 
about  discrimination  coefficients  at  variable  temperatures  and  salinities. 
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4.4.2  Differences  in  Larval  &  Juvenile  Uptake 


We  found  highly  significant  differences  between  uptake  in  larval  and  juvenile  My  a 
arenaria  shells.  In  addition,  larval  and  juvenile  discrimination  coefficients  were  not 
significantly  correlated  for  four  of  the  five  elements  studied,  suggesting  that,  the  phys¬ 
iological  mechanisms  influencing  uptake  are  quite  different  for  the  two  stages.  Not 
only  do  discrimination  coefficients  change  with  ontogeny,  but  the  effects  of  environ¬ 
mental  conditions  on  discrimination  coefficients  change  as  well.  The  one  exception  we 
found  was  for  Ba;  discrimination  coefficients  correlated  positively  between  larvae  and 
juveniles,  suggesting  that  incorporation  of  Ba  into  shell  is  not  as  affected  by  ontogeny 
as  the  other  elements  in  this  study.  This  may  be  because  Ba  is  not  known  to  be  a 
critical  element  in  shell  construction,  nor  is  it  harmful  to  the  organism  when  present 
in  the  extrapallial  fluid.  Ba  concentrations  in  biogenic  carbonate  are  therefore  likely 
to  reflect  ambient  levels. 

The  differences  between  larval  and  juvenile  shell  discrimination  coefficients  are  not 
surprising  given  the  disparate  mechanisms  that  produce  larval  and  juvenile  shell  in 
bivalves.  Production  of  the  first  larval  shell,  or  prodissoconch  I,  starts  by  evagination 
of  a  shell  gland  and  subsequent  spreading  of  the  shell  field.  Prodissoconch  II  formation 
begins  after  the  two  valves  surround  the  body  and  close  against  each  other  (Waller, 
1981).  Early  in  prodissoconch  I  formation,  the  shell  gland  transitions  into  the  mantle 
structure,  which  is  responsible  for  shell  production  thereafter  (Waller,  1981).  These 
intense  morphological  and  developmental  changes  occurring  during  PI  production 
are  likely  to  cause  significant  differences  in  elemental  uptake  compared  to  juvenile 
shell  production  when  the  mantle  has  formed  completely.  Our  results  lend  support 
to  the  idea  that  the  effects  of  physiology  on  larval  and  juvenile  shell  production  are 
sufficiently  different  to  affect  uptake  and  incorporation  of  elements  into  carbonate. 

In  a  study  of  Mg  and  Sr  concentrations  in  nautilus  shells,  Mann  (1992)  found 
consistently  higher  and  more  variable  elemental  concentrations  in  carbonate  formed 
earlier  in  life.  Three  possible  explanations  are  given:  (1)  an  age-related  change  in 
biomineralization,  wherein  elemental  concentrations  depend  on  carbonate  accretion 
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rates;  (2)  maturation  of  the  biomineralization  system  which  results  in  increased  con¬ 
trol  over  shell  chemistry;  or  (3)  larvae  are  exposed  to  stresses  related  to  food  acquisi¬ 
tion  and  protection  that  are  less  problematic  in  older  organisms.  Given  the  evidence 
that  growth  rate  affects  elemental  incorporation  (Onuma  et  ah,  1979;  Purton  et  ah, 
1999;  Sinclair,  2005),  and  the  fact  that  different  processes  are  involved  in  larval  and 
juvenile  shell  accretion,  it  is  plausible  that  any  or  all  of  these  ontogenetic  explanations 
are  applicable  to  our  data. 

Larval  discrimination  coefficients  tended  to  be  closer  to  1  than  juvenile  discrim¬ 
ination  coefficients  for  all  elements  (Fig.  4.4).  This  result  suggests  that  the  more 
advanced  development  of  juveniles  compared  to  larvae  results  in  a  greater  ability  to 
discriminate  against  ions  that  are  not  Ca2+  during  shell  production.  Hirao  et  al. 
(1994)  suggested  that  Pb  content  in  abalone  Haliotis  shells  decreased  due  to  the  de¬ 
velopment  of  a  physiological  mechanism  for  Pb  exclusion  that  becomes  more  efficient 
with  age.  Bivalve  shell  composition  is  influenced  by  the  efficiency  of  the  Ca2+  channel 
that  transports  ions  into  the  extrapallial  fluid  for  shell  construction,  as  well  as  the 
ability  to  discriminate  against  non-Ca2+  ions  (Carre  et  al.,  2006).  Although  there  are 
no  data  published  on  differences  in  channel  function  for  larval  and  juvenile  molluscs, 
large  changes  in  morphology  and  cell  function  during  PI  production  may  influence 
the  composition  of  the  fluid  from  which  the  shell  is  precipitated. 

4.4.3  Conclusions  &;  Implications  for  Future  Studies 

Results  from  studies  such  as  this  one  are  crucial  to  studies  that  rely  on  biogenic  car¬ 
bonate  to  indicate  either  natal  habitat  or  past  environmental  conditions.  Temperature 
and  salinity  effects  on  discrimination  coefficients  of  any  element  can  be  beneficial  for 
connectivity  studies  when  temperature  and  salinity  vary  geographically  while  the  el¬ 
ement  itself  does  not.  However  this  tool  is  useful  only  if  the  study  focuses  on  a  single 
cohort  or  if  the  temperature  and  salinity  conditions  do  not  vary  over  time.  In  con¬ 
trast,  temperature  and  salinity  influences  on  discrimination  coefficients  complicate 
interpretations  when  variation  in  shell  elemental  composition  is  used  as  a  recorder 
of  past  elemental  concentrations  in  water,  or  in  natural  tagging  studies  when  bio- 
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genic  carbonate  formed  outside  of  the  spawning  time  frame  is  used  to  identify  natal 
habitats. 

We  found  that,  temperature  and  salinity  significantly  influenced  shell  chemistry 
of  larval  and  juvenile  My  a  arenaria.  However,  physiological  processes,  whether  in¬ 
fluenced  by  ontogeny  or  environmental  conditions,  have  the  potential  to  complicate 
interpretations  of  elemental  composition  in  biogenic  carbonates.  Variability  in  shell 
chemistry  could  be  interpreted  as  indicating  shifting  environmental  conditions  when 
in  fact  composition  may  reflect  individual  (i.e.  organism-level)  or  ontogenetic  vari¬ 
ability  in  physiological  processes.  Therefore,  care  must  be  taken  to  understand  the 
potentially  confounding  effects  of  physiology  before  shell  chemistry  is  used  in  studies 
as  a  proxy  for  environment  or  to  identify  natal  habitat  for  connectivity  studies. 
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Chapter  5 


Limited  genetic  variation  and 
structure  in  softshell  clams  (My  a 
arenaria)  across  their  native  and 
introduced  range 

5.1  Introduction 

Benthic  marine  habitats  of  the  Northwest  Atlantic  Ocean  (NWA)  are  structured  into 
distinct  biogeographic  provinces  (Engle  and  Summers,  1999).  These  biogeographic 
divisions  are  a  function  of  environmental  gradients  resulting  from  the  synergy  of  the 
coastal  geography  of  Eastern  North  America  with  the  Gulf  Stream  and  Labrador 
Currents,  combined  with  latitudinal  gradients  in  temperature  and  salinity  (Hutchins, 
1947).  The  most  commonly  recognized  biogeographic  divisions  are  the  Nova  Sco¬ 
tian  and  Virginian  Provinces,  with  Cape  Cod  serving  as  the  boundary  between  the 
two  (Hall,  1964;  Hutchins,  1947).  Superimposed  on  these  divisions  is  a  history  of 
Pleistocene  glaciations  that  extirpated  many  benthic  marine  species  from  northern 
latitudes  and  formed  Cape  Cod  (Upham,  1879a, b),  reshaping  regional  patterns  of 
biological  and  genetic  diversity  (Wares  and  Cunningham,  2001;  Wares,  2002;  Hewitt, 
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1996). 


The  presence  of  distinct  biogeographic  provinces  in  the  NWA  has  significant  im¬ 
plications  for  management  of  fish  and  invertebrates  in  this  region  because  species 
spanning  multiple  provinces  of  the  NWA  can  have  populations  adapted  to  local  envi¬ 
ronmental  conditions.  For  example,  the  Atlantic  Silverside  Menidia  menidia  exhibits 
heritable  local  variation  in  growth  rate  and  vertebral  number,  resulting  in  a  latitudinal 
phenotypic  cline  across  the  NWA  (Present  and  Conover,  1992;  Billerbeck  et  ah,  1997; 
Yamahira  et  ah,  2006).  On  a  smaller  scale,  the  mussel  Mytilus  edulis  exhibits  a  sharp 
cline  in  the  leucine  aminopeptidase  (LAP)  allele  across  salinity  gradients  in  Long 
Island  Sound  (Gardner  and  Kathiravetpillai,  1997;  Gardner  and  Palmer,  1998).  The 
presence  of  regional  genetic  structure,  particularly  if  it  is  locally  adaptive,  needs  to 
be  accounted  for  in  fisheries  management  so  that  genetic  diversity  is  conserved  and 
locally  adaptive  gene  complexes  are  not  disrupted  through  indiscriminate  stocking 
(Hansen,  2002). 

My  a  arenaria  is  a  commercially  important  bivalve  with  a  contemporary  distribu¬ 
tion  that  includes  1)  the  northwest  Atlantic  ranging  from  Nova  Scotia  to  Virginia,  2) 
the  North  Sea  and  European  waters,  including  the  Black,  Baltic,  Wadden,  White,  and 
Mediterranean  Seas,  and  3)  northeast  Pacific  from  San  Francisco  to  Alaska  (Strasser, 
1999).  M.  arenaria  has  a  complex  history  of  extensive  global  distributions,  with 
several  extinctions  and  re-colonization  events  (reviewed  in  Strasser  1999).  M.  are¬ 
naria  originated  in  the  Pacific  Ocean  during  the  Miocene  then  extended  its  range  to 
the  Atlantic  and  European  waters  in  the  early  Pliocene.  Extinction  of  Pacific  and 
European  populations  in  the  early  Pleistocene  left  the  only  surviving  populations  in 
the  NWA  until  recent  history  (MacNeal,  1965).  M.  arenaria  re-invaded  European 
waters  in  the  17th  century  after  being  brought  from  the  NWA  by  Vikings  (Petersen 
et  ah,  1992).  In  the  late  19th  century  M.  arenaria  was  reintroduced  into  the  Pacific, 
first  accidentally  then  as  a  potential  commercial  fishery  (Carlton,  1979;  Powers  et  ah, 
2006).  The  natural  and  introduced  distribution  of  M.  arenaria  results  partly  from  the 
clam’s  ability  to  withstand  wide  salinity  and  temperature  ranges,  and  its  capability 
of  inhabiting  different  sediment  types  from  fine  mud  to  coarse  sand  (Newell  and  Hidu, 
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1982;  Abraham  and  Dillon,  1986;  Hidu  and  Newell,  1989). 

Despite  characteristics  that  made  softshell  clams  a  successful  invasive  species  in 
European  waters  and  the  northeast  Pacific,  the  last,  two  decades  have  seen  appreciable 
declines  in  softshell  clam  landings  in  New  England  (Brousseau.  2005;  Anonymous, 
2007).  This  decline  has  been  attributed  to  habitat  degradation  or  loss,  overfishing, 
contamination,  and  predation  by  invasive  species  (Brousseau,  2005).  Managers  and 
state  agencies  have  enacted  various  management  strategies  to  combat  these  declines, 
including  using  protective  nets  to  reduce  predation  on  newly  recruited  clams,  and 
seeding  flats  using  hatchery-reared  juveniles  (Marcotti  and  Leavitt,  1997,  H.  Lind, 
pers.  comm.). 

While  stocking  of  fish  and  shellfish  is  a  long-standing  pract  ice,  research  is  increas¬ 
ingly  showing  that  the  genetic  impacts  of  stocking  cannot  be  ignored.  Stocking  should 
seek  to  maintain  levels  of  genetic  diversity  (Waples  and  Do,  1994);  although  multiple 
individuals  are  spawned  to  produce  seed  clams,  it  is  unknown  whether  the  genetic 
diversity  represented  among  these  individuals  is  reduced  in  comparison  to  naturally 
occurring  cohorts,  where  entire  adult,  populations  spawn  simultaneously  (Brousseau, 
1978).  In  addition,  given  that,  brood  stock  is  not,  always  taken  from  the  flat,  into 
which  seed  clams  are  stocked,  locally  appropriate  genotypes  could  be  introduced  into 
inappropriate  areas.  For  example,  Mya  arenama  exhibit,  local  variation  in  resistance 
to  paralytic  shellfish  toxins  (Connell  et.  al.,  2007).  Seeding  flats  using  brood  stock 
from  other  clam  populations  may  result  in  either  reduction  of  the  locally  dominant, 
alleles  due  to  success  of  the  introduced  seed  clams,  or  significant  loss  of  seed  clams 
due  to  a  lack  of  a  genetic  background  appropriate  to  the  local  environment.  Similar 
declines  in  local  fitness  have  been  documented  in  salmonids  (Hansen,  2002). 

Previous  genetic  studies  on  softshell  clams  have  found  limited  genetic  diversity 
despite  the  wide  geographic  ranges  represented  among  studies.  Morgan  et  al.  (1978) 
used  allozymes  to  study  Mya  arenaria  genetic  variation  in  the  NWA,  and  found  low 
polymorphism  and  low  heterozygosity  per  individual  for  both  populations  examined. 
Caporale  et.  al.  (1997)  found  similar  low  variability  in  three  regions  of  the  NWA  (12 
locations  total)  using  the  internal  transcribed  spacer  ribosomal  DNA  region  (nDNA), 
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and  concluded  that  although  the  three  regions  were  not  genetically  distinct,  the  data 
from  the  study  were  insufficient  to  indicate  a  panmictic  population.  More  recently, 
Lasota  et  al.  (2004)  used  allozymes  to  study  seven  locations  in  the  northeast  Atlantic 
and  two  in  the  North  Sea.  They  also  found  low  genetic  variability  and  a  lack  of  genetic 
differentiation,  and  concluded  that  M.  arenaria  is  a  successful  invader  despite  a  high 
degree  of  genetic  homogeneity.  They  suggested  the  patterns  observed  were  evidence 
of  rapid  population  expansion,  allele  neutrality,  and  high  gene  flow.  However,  nDNA 
is  known  to  evolve  slower  than  mitochondrial  DNA  (mt.DNA),  and  allozyme  studies 
may  mask  underlying  sequence  variation.  Therefore  the  results  seen  in  these  studies 
might  be  because  of  the  markers  chosen  to  conduct  the  studies. 

In  this  study,  we  examine  population  genetic  variability  of  Mya  arenaria  across  its 
natural  range  in  the  NWA  and  portions  of  its  introduced  range  in  the  northeast  Pa¬ 
cific  and  European  waters  using  the  highly  variable  mitochondrial  cytochrome  oxidase 
I  (COI)  gene  that  commonly  resolves  phylogeographic  structure  in  marine  inverte¬ 
brates  (Wares,  2002;  Barber  et  ah,  2006)  including  bivalves  (King  et  ah,  1999;  May 
et  ah,  2006).  First,  we  examine  how  populations  may  be  geographically  structured 
across  the  NWA  to  determine  whether  the  distinct  environments  and  biogeographic 
provinces  partition  softshell  clams  into  genetically  distinct  regional  stocks.  Second, 
we  compare  NWA  to  NEP  and  NSE  populations  to  examine  the  geographic  origins 
of  these  populations  and  the  effects  of  recent  introduction  on  genetic  diversity.  The 
results  of  this  study  have  implications  for  management  of  softshell  clams  in  New  Eng¬ 
land,  in  addition  to  insights  gained  about,  historical  extinction  and  colonization  events 
of  M.  arenaria  with  reference  to  biogeographic  boundaries  and  glaciation. 

5.2  Methods 

Juvenile  and  adult  M.  arenaria(N  =  212)  were  collected  between  2001  and  2006  from 
12  locations:  one  northeast  Pacific  site  ( n  =  20),  ten  NWA  sites  (n  =  177),  and  one 
North  Sea,  Europe  site  (n  =  15)  (Fig.  5.1,  Table  5.1).  Most  M.  arenaria  were 
frozen  after  collection  to  prevent  DNA  degradation,  and  then  transferred  to  70-95% 
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Figure  5.1:  Unrooted  minimum-spanning  tree  depicting  the  relationship  of  the  27 
mitochondrial  COI  haplotypes  from  217  Mya  arenana  individuals,  collected  from  12 
sites  in  the  NWA  (white;  n  =  177),  one  site  in  the  North  Sea,  Europe  (black;  n  —  15), 
and  one  site  in  the  northeast  Pacific  (gray;  n  =  20).  Line  distance  between  circles 
corresponds  to  the  number  of  nucleotide  differences  (1  or  2);  size  and  number  of 
divisions  in  each  haplotype  circle  correspond  to  the  number  of  individuals. 
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ethanol  for  at.  least  24  hours  prior  to  DNA  extraction  to  improve  the  success  of  DNA 
extractions.  Some  individuals  were  preserved  directly  in  ethanol  without  freezing. 
For  clams  <  1cm  total  length,  we  used  the  entire  clam  for  DNA  extraction.  For 
larger  clams  we  extracted  DNA  from  small  fibers  of  adductor  muscle  tissue.  All  DNA 
extractions  were  performed  with  a  10%  Chelex®(BioRad)  solution  (Walsh  et  al., 
1991).  A  661  bp  fragment  of  the  mitochondrial  cytochrome  oxidase  subunit-I  gene 
(COI)  was  amplified  via  polymerase  chain  reaction  (PCR)  using  the  primers  HCO- 
2198  and  LCO-1490  (Folmer  et  ah,  1994).  PCR  occurred  in  25  1  reactions  with  2.5  1 
of  lOx  buffer,  2  1  MgCb  (25  mM),  2.5  1  DNTPs  (8  mM),  1.25  1  of  each  10  mM  primer, 
1  1  of  template,  and  0.625  units  of  Amplitaq  (Perkin  Elmer).  Hot-start  thermocycling 
parameters  were  as  follows:  initial  denaturation  94°C  (3  min);  followed  by  38  cycles 
of  94°C  (30  s),  50°C  (30  s),  72°C  (45  s);  then  a  final  extension  of  72°C  (10  min). 

PCR,  products  were  visualized  on  1%  agarose  PAC  1%  sodium  hydroxide  and  boric- 
acid  gels,  and  then  enzymatically  prepared  for  sequencing  by  digestion  in  0.5  units 
of  Shrimp  Alkaline  Phosphotase  and  5  units  of  exonuclease  per  5  1  of  PCR  product, 
incubated  at  37°C  for  30  min  followed  by  80°C  for  15  min.  PCR  product  was  cleaned 
using  isopropanol  precipitation,  and  sequencing  reactions  were  performed  for  both  for¬ 
ward  and  reverse  strands  on  an  ABI  377  (Applied  Biosystems)  using  BigDye(c)( Perkin 
Elmer)  terminator  chemistry.  Complementary  strands  for  each  sample  were  proof¬ 
read  and  aligned  in  Sequencher7  M,  and  translations  confirmed  using  MaClade  4.05 
(Maddison  and  Maddison,  2002)  . 

To  explore  regional  distribution  of  genetic  diversity  in  Mya  arenaria ,  we  calculated 
haplotype  diversity  (/i),  nucleotide  diversity  (7r),  and  theta  (6)  for  all  populations 
using  Arlequin  3.1  (Excoffier  et  ah,  2005).  To  explore  patterns  of  phylogeographic 
structure  we  constructed  a  minimum-spanning  tree  using  the  MINSPNET  algorithm 
as  employed  in  Arlequin.  Frequency  of  haplotypes  was  then  plotted  against  geography 
for  NWA  populations. 

To  further  explore  geographic  genetic  structure,  we  investigated  genetic  partitions 
in  analysis  of  molecular  variance  (AMOVA)  as  implemented  in  Arlequin.  Values  of 
4>st  were  calculated  with  statistical  significance  determined  by  20.000  random  per- 
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mutations.  Analyses  were  run  both  unstructured,  and  structured  into  two  regions 
(NWA  +  NEP  and  NSE)  or  three  regions  (NWA,  NEP,  NSE).  Patterns  of  genetic 
structure  were  similarly  estimated  within  the  NWA  by  excluding  NEP  and  NSE  pop¬ 
ulations.  NWA  analyses  were  both  unstructured  and  assuming  three  regions:  north 
of  Cape  Cod,  Cape  Cod,  and  South  of  Cape  Cod.  To  further  examine  patterns  of 
genetic  exchange,  pairwise  F st  values  were  calculated  among  all  populations  with 
20,000  permutations  used  to  establish  significance. 


Because  of  extremely  low  levels  of  genetic  variation  within  the  dataset,  we  tested 
for  neutrality  by  calculating  Fu’s  F5  statistics  (Fu,  1997),  which  establishes  whether 
non-neutrality  might  be  due  to  population  growth  and  range  expansion.  To  further 
explore  the  possibility  of  recent  demographic  or  spatial  population  expansion,  we 
used  mismatch  distributions,  which  compares  the  expected  and  observed  number  of 
differences  between  pairs  of  haplotypes  (Rogers  and  Harpending,  1992;  Ray  et  al., 
2003).  Finally,  we  used  Bayesian  Markov  Chain  Monte  Carlo  analysis  of  molecular 
sequences  to  produce  a  Bayesian  skyline  plot  using  BEAST  vl.4  and  Tracer  vl.4 
(Drummond  and  Rambaut,  2006),  which  plots  population  size  over  time  and  estimates 
the  approximate  time  since  population  expansion  (Drummond  et  ah,  2002,  2005).  We 
used  MODELTEST  (Posada  and  Crandall,  1998)  implemented  in  PAUP*  ver.4.0bl0 
(Swofford,  1998)  to  find  the  most  appropriate  model  for  BEAST  (Hasegawa,  Kishino, 
and  Yano  Model).  A  strict  molecular  clock  was  used  to  produce  the  skyline  plot, 
which  was  based  on  five  skyline  groups.  We  ran  the  program  using  default  priors  for 
Bayesian  skyline  analysis  for  50  million  generations,  and  repeated  the  program  run 
four  times  to  increase  effective  sample  size  and  assure  that  results  were  converging. 
Results  reported  in  mutational  units  were  converted  to  years  for  the  skyline  plot  by 
assuming  a  molluscan-specific  COI  divergence  rate  of  either  1%  per  million  years 
(%Myr-1)  for  all  COI  sites  or  5  %Myr-1  for  third  positions  alone  (Marko,  2002). 
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5.3  Results 


A  total  of  661  bp  of  COI  was  collected  from  212  individuals,  yielding  only  27  unique 
haplotypes  that  differed  by  one  or  two  nucleotide  substitutions,  all  in  the  third  codon 
position  and  silent  with  one  exception,  haplotype  I  (Table  5.2).  There  was  one  dom¬ 
inant  haplotype  (A)  found  at  all  of  the  locations  sampled,  ranging  in  frequency  from 
0.65  to  1.00  for  individual  populations,  with  an  overall  frequency  of  0.79.  Of  the 
remaining  26  haplotypes,  only  five  were  found  more  than  once  in  a  single  population, 
ranging  in  frequency  from  0.10  to  0.27  (haplotypes  B-F).  Two  private  alleles,  haplo¬ 
types  that  occur  more  than  one  time  in  only  one  site  (Slatkin,  1985),  were  found  in 
NSE  (Haplotype  E)  and  NEP  (Haplotype  G). 

Haplotype  diversity  (/?.)  in  the  NWA  ranged  from  0.178  to  0.648  (Table  5.1).  Com¬ 
parable  levels  of  haplotype  diversity  occurred  in  NSE  ( h  =  0.65)  and  NEP  (/?  =  0.57) 
populations.  Nucleotide  diversity  was  low  for  all  NWA  populations,  ranging  from 
0.0003  to  0.001  (Table  5.1),  while  tt  =  0.0010  in  NSE  and  ir  =  0.0012  in  NEP.  Theta 
ranged  from  1.854  to  0.549  in  the  NWA  and  was  0.615  and  1.127  in  NSE  and  NEP, 
respectively  (Table  5.1).  There  were  no  clear  geographic  patterns  in  genetic  diversity 
measures. 

Consistent  with  the  low  nucleotide  diversity,  the  minimum  spanning  tree  of  M. 
arenaria  COI  haplotypes  revealed  a  star-shaped  phylogeny  (Fig.  5.1).  The  dominant 
haplotype  (A)  was  located  at  the  center  of  the  star  with  21  of  26  remaining  haplo¬ 
types  differing  from  haplotype  A  by  a  single  nucleotide  substitution.  Five  haplotypes 
differed  by  2  mutational  steps  (haplotypes  J,  O.  T,  U.  AA).  No  geographic  structure 
is  evident  in  the  minimum  spanning  tree  topology  and  NEP  and  NSE  haplotypes  are 
scattered  throughout  the  tree.  Plotting  the  frequency  of  the  6  non-singleton  haplo¬ 
types  revealed  no  clear  phylogeographic  patterns  in  the  NWA  except  for  the  lack  of 
genetic  diversity  in  Nova  Scotia  (Fig.  5.2). 

Results  from  AMOVA  found  the  majority  of  variability  was  within  populations,  re¬ 
gardless  of  any  structure  imposed  on  the  locations  sampled  (Table  5.3).  Unstructured 
AMOVA  analyses  indicated  the  presence  of  subtle  but  significant  genetic  structure 
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_ Location _ 

Haplotvpe  NS  NB  ME  OMA  BMA  MM  A  WMA  NY  MMD  EMD  NEP  NSE  ALL 

A  1.00  (20)  0.73  (8)  0.86  (18)  0.91(20)  0.9  (17)  0.75  (15)  0.76  (19)  0.73  (11)  0.80  (12)  0.67  (6)  0.65  (13)  0.53  (8)  0.788  (167) 

B  0.09  (1)  0.04  (1)  0.07  (1)  0.07  (1)  0.05  (1)  0.20  (3)  0.038  (8) 
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Table  5.2:  Haplotype  distributions  for  Mya  arenaria.  Haplotype  frequencies  are  given  for  each  locality  sampled,  with  the  number 
of  individuals  per  haplotype  in  parentheses.  ALL  column  is  sum  of  all  sites  sampled. 
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Figure  5.2:  Distribution  of  mitochondrial  COI  haplotypes  for  Mya  arenana  in  the 
NWA.  Gray  shades  are  unique  haplotypes  found  in  only  one  location:  patterns  are 
haplotypes  shared  among  two  or  more  locations.  See  Table  5.1  for  site  abbreviations 
and  sample  sizes. 
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Among  regions 

Among  populations  within  regions 

Within  populations 

df 

Var. 

%  Var. 

4Vt 

df 

Var. 

%  Var. 

4>sc 

df 

Var. 

%  Var. 

Osi 

All  Locations 

Unstructured 

11 

0.0062 

2.67 

200 

0.2271 

97.3 

0.0267* 

Structured:  2  groups 

1 

0.0429 

15.9 

0.159 

10 

0.00003 

0.01 

0.0001 

200 

0.2271 

84.1 

0.1591* 

Structured:  3  groups 

2 

0.0239 

9.57 

0.957 

9 

-0.0014 

-0.54 

-0.0060 

200 

0.2271 

91.0 

0.0903* 

NWA  Only 

Unstructured 

9 

0.00008 

0.04 

167 

0.2020 

99.9 

0.0004 

Structured 

2 

0.00043 

0.21 

0.0021 

7 

-0.0002 

-0.12 

-0.0012 

167 

0.2020 

99.9 

0.0010 

Table  5.3:  Results  of  AMOVAs  to  determine  the  source  of  genetic  variation  in  Mya 
arenaria  COI.  *  indicates  that  the  value  is  significant  at  the  p  =  0.005  level.  See  text 
for  group  descriptions. 

( 4>st  =  0.0267,  p  <  0.005)  with  3%  of  the  variation  between  populations  and  97% 
of  the  variation  within  populations.  Imposing  regional  partitions  comparing  North 
American  (NWA+NEP)  and  European  (NSE)  populations  produced  (psr  =  0.159 
(p  <  0.005)  with  15.9%  of  variation  among  regions,  0.01%  among  populations  within 
regions,  and  84.1%  of  the  variation  within  populations.  When  we  imposed  three  re¬ 
gional  partitions,  NWA,  NEP,  and  NSE,  (psr  =  0.0903  ( p  <  0.005)  with  91%  of  the 
variation  within  populations,  9%  among  regions,  and  no  variation  among  populations 
within  regions. 

Unstructured  AM  OVA  analyses  across  the  NWA  revealed  no  significant  genetic 
structure  with  99.9%  of  all  genetic  variation  contained  within  populations  (Table  5.3). 
Similarly,  when  locations  were  grouped  into  regions  north  of  Cape  Cod,  Cape  Cod, 
and  south  of  Cape  Cod  (Table  5.1),  (j)sr  =  0.0010  (n.s.)  with  99.9%  of  all  genetic 
variation  contained  within  populations  and  no  significant  variation  among  regions  or 
among  populations  within  regions. 

Population  pairwise  Fst’s  were  significant  and  high  between  NSE  and  7  popula¬ 
tions,  including  NEP  (FgT-  =  0.09  to  0.222;  p  <  0.05)  (Table  5.4).  In  addition,  NEP 
showed  moderate  divergence  with  QMA  (F st  =  0.039;  p  =  0.006).  Within  the  NWA, 
sites  showed  low  variability  between  populations.  Only  the  northernmost  population 
(Nova  Scotia)  had  significant  F st  values,  with  moderate  divergence  between  Nova 
Scotia  and  populations  from  New  Brunswick  (Fst  =  0.058)  and  Eastern  Bay  Mary¬ 
land  (Fsx  =  0.097).  In  addition,  there  was  minimal  divergence  (F5r  =  0.02)  between 
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NSE 

NS 

NB 

ME 

OMA 

BMA 

MMA 

WMA 

NY 

MMD 

EMD 

NEP 

NSE 

- 

0.0007 

0 

0.0003 

0.001 

0.0003 

0.009 

0.008 

NS 

0.222* 

- 

0.036 

0.026 

0.023 

NB 

0.065 

0.058* 

ME 

0.143* 

-0.002 

0.008 

- 

QMA 

0.182* 

-0.004 

0.028 

0.001 

- 

0.006 

BMA 

0.166* 

0.003 

0.006 

-0.001 

0.001 

- 

MMA 

0.115* 

0.000 

-0.018 

0.001 

0.003 

-0.008 

- 

WMA 

0.099* 

0.000 

-0.010 

-0.010 

0.005 

0.001 

-0.001 

- 

NY 

0.084 

0.02* 

-0.023 

0.003 

0.012 

0.007 

-0.013 

-0.015 

- 

MMD 

0.092 

0.020 

-0.023 

0.000 

0.009 

0.005 

-0.003 

-0.010 

-0.019 

EMD 

0.052 

0.097* 

-0.024 

0.012 

0.060 

0.047 

0.013 

-0.020 

-0.010 

-0.006 

NEP 

0.104* 

0.045 

-0.023 

0.013 

0.039* 

0.022 

0.009 

-0.001 

0.013 

0.012 

-0.015 

Table  5.4:  Pairwise  population  comparisons,  F St  (below  diagonal)  and  their  associ¬ 
ated  p  values  if  significant  (above  diagonal).  *  indicates  that  the  value  is  significant 
at  the  p  —  0.05  level. 

Nova  Scotia  and  New  York.  However,  no  other  pairwise  Fsr’s  were  significant  at  the 
p  <  0.05  level. 

Some  haplotypes  found  in  NSE  and  NEP  were  not  shared  with  the  NWA.  To 
determine  whether  this  was  due  to  inadequate  sampling,  we  constructed  a  rarefaction 
curve  for  the  NWA  (Fig.  5.3)  using  equations  appropriate  for  population  sample  sizes 
much  smaller  than  total  sample  size  (Heck  et  al.,  1975).  The  rarefaction  curve  did  not 
reach  an  asymptote  over  the  range  of  the  number  of  individuals  sampled.  Although 
it  is  not  valid  to  extrapolate  the  curve  and  predict  how  many  haplotypes  might  be 
present  in  the  NWA,  the  slope  of  the  curve  does  appear  to  be  leveling  off,  which 
suggests  that  a  majority  of  the  available  haplotypes  were  sampled. 

Fu’s  F5  Statistic  was  significantly  large  and  negative  for  9  of  the  12  populations 
(Table  5.1)  suggesting  non-equilibrium  dynamics.  Mismatch  analysis  revealed  no 
significant  fit  to  the  models  of  recent  spatial  or  demographic  expansion  (data  not 
shown).  However,  p- values  were  generally  lower  (i.e.  the  model  was  a  better  fit)  for 
the  spatial  expansion  hypothesis. 

Support  for  a  range  expansion  comes  from  the  Bayesian  skyline  plot  indicating 
that  Mya  arenaria  populations  in  the  NWA  were  much  smaller  in  recent  history 
(Fig.  5.4).  The  plot  indicates  that  a  pronounced  demographic  expansion  event  took 
place  in  NWA  populations  of  M.  arenaria  between  75,000  and  15,000  years  ago.  These 
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Figure  5.3:  Rarefaction  curve  constructed  using  data  from  Mya  arenaria  populations 
in  the  NWA. 

values  correspond  to  0.00035  mutational  units  and  a  mutation  rate  of  0.005  to  0.025 
mutations  Myr-1,  based  on  the  clock  calibrations  of  Marko  (2002). 


5.4  Discussion 

5.4.1  Patterns  in  the  Northwest  Atlantic 

Genetic  analysis  of  Mya  arenaria  populations  across  the  Northwest  Atlantic  revealed 
a  near  complete  absence  of  genetic  structure.  This  result  stands  in  contrast  to  previ¬ 
ous  studies  of  other  marine  species  that  show  pronounced  phylogeographic  structure 
in  the  NWA  (see  Wares  2002  for  a  review),  particularly  among  populations  along 
the  northern  and  southern  coastline  of  the  NWA  (e.g.  Waldman  et  al.,  1996;  Smith 
et  ah,  1998;  Dahlgren  et  ah,  2000;  Brown  et  ah,  2001).  Within  the  NWA  the  only 
evidence  for  differentiation  among  northern  and  southern  populations  comes  from 
three  significant  pairwise  F st  values  among  Nova  Scotia  populations  and  those  to  the 
south.  These  significant  differences  are  likely  driven  by  the  lack  of  diversity  in  the 
NS  sample,  where  only  the  dominant  haplot.ype  was  observed. 

The  lack  of  genetic  diversity  and  limited  genetic  structure  reported  here  echoes 
previous  genetic  studies  on  softshell  clams  (Morgan  et  ah,  1978;  Caporale  et  ah,  1997; 
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Figure  5.4:  Bayesian  skyline  plot  derived  from  Mya  arenaria  NWA  sequences.  The 
solid  line  is  the  median  estimate  of  population  size,  and  the  shaded  region  shows 
95%  highest  posterior  density  limits  (see  Drummond  et.  ah,  2005).  The  dashed  line 
indicates  where  in  time  the  population  expanded.  Axis  A  is  the  years  before  present 
when  a  1%  per  million  years  divergence  rate  is  used;  Axis  B  is  the  years  before  present 
when  a  5%  per  million  years  divergence  rate  is  used  (see  text  for  details).  The  gray 
arrows  on  each  axis  represent  the  approximate  timing  of  the  last  glacial  maximum. 
Population  size  on  the  y-axis  is  relative  to  its  current  size. 
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Lasota  et  al.,  2004).  The  concordant  results  among  these  multiple  studies  provide 
strong  evidence  for  lack  of  genetic  boundaries  in  Mya  arenaria.  The  lack  of  genetic 
structure  among  NWA  populations  of  M.  arenaria  could  result  from  high  levels  of 
gene  flow,  as  suggested  by  Lasota  et  al.  (2004),  due  to  extensive  dispersal  and  open 
population  dynamics.  This  species  has  a  planktonic  larval  phase  that  can  last  up  to 
three  weeks  in  the  water  column,  during  which  time  the  larva  feeds  on  algae  and  is 
transported  by  currents  (Abraham  and  Dillon,  1986).  Transport  via  strong  currents 
along  the  NWA  could  promote  high  dispersal  and  gene  flow  among  NWA  populations. 
Genetic  mixing  could  be  further  augmented  by  human-mediated  transport  within  the 
NWA.  High  gene  flow,  however,  should  increase  the  effective  population  size,  reducing 
drift  and  resulting  in  high  levels  of  genetic  diversity.  This  expectation  is  contradicted 
by  the  minimum  spanning  tree  and  the  minimal  genetic  diversity  measures.  Thus, 
while  high  dispersal  may  contribute  to  genetic  homogeneity  in  M.  arenaria ,  other 
processes  must  also  be  acting  to  reduce  diversity  within  populations. 

One  process  that  could  contribute  to  limited  genetic  structure  across  the  NWA 
is  a  recent  population  expansion  event,  also  suggested  by  Lasota  et  al.  (2004).  The 
NWA  was  heavily  impacted  by  ice  sheets  during  the  Pleistocene  glacial  period  ap¬ 
proximately  2.5  million  years  ago,  with  the  southern  limit  of  glaciation  near  Cape 
Cod,  Massachusetts  (Shackleton  et  al.,  1984;  Cronin,  1988).  It  is  hypothesized  that 
the  ice  sheet  caused  a  southward  shift  in  population  ranges  of  NWA  intertidal  species. 
After  glaciers  subsided  approximately  18,000  years  ago,  individuals  from  the  south 
spread  into  previously  unavailable  northern  habitats  (Wares  and  Cunningham,  2001; 
Wares,  2002).  These  re-invasions  may  manifest  themselves  in  genetic  structure  that 
suggests  population  expansion,  i.e.  limited  genetic  structure  across  areas  affected  by 
glaciation.  Wares  (2002)  compiled  a  comprehensive  map  of  genetic  variability  in  NWA 
intertidal  species  and  found  substantial  evidence  of  recolonization  of  ice-covered  areas. 
Generally  modern  populations  located  north  of  the  southern  extent  of  the  ice  sheet 
have  lower  genetic  diversity  than  those  located  to  the  south,  with  shifts  in  genetic 
structure  located  near  Cape  Cod. 

Given  the  contemporary  distribution  of  Mya  arenaria  and  the  geographic  extent 
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of  Pleistocene  glaciations,  northern  populations  in  the  NWA  likely  represent  a  recent 
range  expansion.  Evidence  for  a  demographic  expansion  comes  from  the  star-like  phy- 
logeny,  low  genetic  diversity  measures,  and  the  significantly  large  and  negative  values 
of  Fu’s  Fs.  Further  evidence  of  demographic  expansion  comes  from  the  Bayesian 
skyline  plot  produced  using  NWA  data  (Fig.  5.4).  The  plot  shows  that  M.  arenaria 
populations  in  the  NWA  expanded  between  15,000  and  75,000  years  ago  depending 
on  mutation  rate.  Pleistocene  ice  sheets  retreated  beginning  approximately  18,000 
years  ago,  which  is  within  this  range.  Although  the  data  from  Bayesian  Skyline  plot 
analyses  do  not  produce  exact  estimates  without  stringent  assumptions  regarding 
mutation  rates  and  molecular  clocks,  the  data  do  not  contradict  the  hypothesis  of 
post-glacial  expansion. 


Given  a  northward  range  expansion  following  glacial  retreat,  one  might  predict 
decreasing  genetic  diversity  measures  with  increasing  latitude.  Although  the  north¬ 
ernmost  population  in  Nova  Scotia  had  the  lowest  diversity  measures,  there  was  no 
pattern  across  the  remainder  of  the  NWA.  This  result  suggests  that  while  dispersal 
and  gene  flow  may  be  high  among  most  populations  spanning  the  NWA,  Nova  Scotia 
populations  may  be  relatively  isolated. 


An  alternative  explanation  for  the  signal  of  demographic  expansion  is  recovery 
following  a  selective  sweep,  where  selectively  advantageous  haplotypes  go  to  fixation 
(e.g.  Berry  et.  al.,  1991).  Recovery  from  a  recent  selective  sweep  could  also  yield  a  star- 
like  phylogeny  and  lower  genetic  diversity  that  would  inflate  genetic  similarity  and 
gene  flow  estimates  among  populations.  While  a  selective  sweep  cannot  be  excluded, 
given  the  signal  of  demographic  expansion  during  the  Pleistocene  when  demographic 
expansion  would  be  required  to  achieve  contemporary  distributions,  a  selective  sweep 
may  be  a  less  parsimonious  explanation. 
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5.4.2  Patterns  across  the  Northwest  Atlantic,  European  wa¬ 
ters  and  the  northeast  Pacific 

While  populations  in  the  NWA  had  minimal  genetic  structure,  the  strongest  signal  of 
regional  genetic  structure  comes  from  comparing  NWA  populations  to  NSE.  Struc¬ 
tured  AMOVA  results  with  NSE,  NEP,  and  NWA  defined  as  separate  regions  resulted 
in  a  significant  <$>st  of  0.0903  (Table  5.3).  Furthermore,  of  11  pair-wise  comparisons,  a 
total  of  7  pair-wise  F$t  values  among  NSE  and  the  NWA  were  significant  (Table  5.2). 
This  result  indicates  that  despite  being  introduced  from  NWA  populations,  there  are 
significant  genetic  differences  among  these  regions. 

Surprisingly,  NEP  and  NSE  populations  also  exhibited  the  highest  levels  of  ge¬ 
netic  diversity.  Previous  studies  have  found  higher  than  expected  genetic  diversity 
in  recently  colonized  habitats  due  to  multiple  invasions  that  reduced  the  founder  ef¬ 
fect  (e.g.  Stepien  et  ah,  2002;  Voisin  et  ah,  2005;  Simon-Bouhet.  et  ah,  2006).  In 
these  studies,  however,  the  diversity  observed  was  not  higher  than  that  of  the  re¬ 
gion  of  origin,  as  was  the  case  in  our  study.  The  unusually  high  levels  of  haplotypes 
and  nucleotide  diversity  observed  in  the  introduced  North  Sea  population  (Table  5.1) 
potentially  suggests  introductions  from  multiple  source  populations,  a  result  that  is 
consistent  with  significant  F st  values  among  NSE  and  some  NWA  populations  but 
no  others  (Table  5.3). 

Curiously,  both  Pacific  and  European  populations  also  contained  private  haplo¬ 
types  (Fig.  5.1).  The  presence  of  unique  haplotypes  found  multiple  times  in  a  single 
population  suggests  genetic  isolation  (Hartl  and  Clark,  1997).  Given  the  geographic 
separation  of  the  NWA,  Pacific  and  European  waters,  observation  of  genetic  isolation 
should  be  expected.  However,  this  result  is  surprising  given  that  both  NEP  and  NSE 
populations  are  thought  to  have  been  introduced  from  the  NWA  within  the  last  150 
and  400  years,  respectively.  This  seems  a  particularly  short  amount  of  time  for  local 
variation  to  both  evolve  in  situ  and  increase  in  frequency  sufficiently  to  be  detected 
by  sampling  15-20  individuals.  In  contrast,  no  private  haplotypes  were  detected  in 
sampling  of  177  individuals  from  the  entire  range  of  M.  arenaria  in  the  NWA. 
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As  a  heuristic  we  constructed  a  rarefaction  curve  plotting  number  of  haplotypes 
versus  number  of  samples.  Although  the  slope  shallows,  it  did  not  asymptote  over 
the  range  of  number  of  individuals  sampled  (Fig.  5.3)  indicating  that,  sampling  177 
individuals  was  insufficient  in  the  NWA  to  detect  the  rare  haplotypes  that  founded 
NEP  and  NSE  populations.  Paradoxically,  if  they  are  in  very  low  frequencies  in  the 
NWA,  it  seems  unlikely  that  they  would  be  introduced  to  the  Pacific  and  European 
waters.  One  interpretation  of  this  result  is  that  these  private  haplotypes  may  represent 
ancestral  polymorphism  from  relic  populations  that  survived  the  extinction  events  in 
the  Pacific  and  European  waters.  However  if  these  were  relic  haplotypes,  genetic 
divergence  in  excess  of  one  mutational  step  would  be  expected,  suggesting  that  these 
are  indeed  introductions  of  rare  NWA  haplotypes. 

5.4.3  Management  Implications 

One  of  the  current  management  strategies  for  NWA  softshell  clam  populations  is  to 
increase  local  abundances  by  seeding  flats  with  hatchery-reared  juvenile  clams.  As 
has  been  demonstrated  in  fish,  this  approach  has  the  potential  to  decrease  or  alter 
genetic  variability  by  introducing  non-native  genotypes  that  may  affect  the  fitness  of 
both  introduced  and  native  stocks  (Hansen,  2002). 

The  low  genetic  diversity  and  minimal  genetic  structure  observed  in  COI  combined 
with  previous  results  showing  limited  genetic  diversity  in  Mya  arenaria  using  nuclear 
sequences  (Caporale  et,  al.,  1997)  and  allozvmes  (Morgan  et.  al.,  1978;  Lasota.  et  ah, 
2004)  suggests  that  the  brood  stock  origins  may  not  be  critical  to  maintaining  levels 
of  genetic  diversity  and  patterns  of  genetic  structure  across  the  NWA.  Results  of  t  his 
study  suggest  that  brood  stocks  should  be  quite  similar  regardless  of  their  locality, 
and  their  resulting  juvenile  seed  clams  are  likely  interchangeable  across  geography. 

Although  we  did  not  detect  genetic  structure  using  the  mitochondrial  COI  gene, 
there  may  yet  be  other  genes  that  might  show  variability  within  the  NWA.  However, 
given  that  Mya  arenaria  is  recovering  from  a  severe  reduction  in  genetic  diversity,  the 
odds  of  detecting  neutral  genetic  variation  that  correspond  to  locally  adaptive  gene 
complexes  is  remote.  Given  that  local  adaptation  has  been  noted  in  the  softshell  clam 
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for  toxin  resistance  (Connell  et  al.,  2007),  there  may  very  well  be  important  regional 
genetic  differences  among  clam  stocks  in  non-neutral  genes.  As  such,  seeding  from 
local  stocks  should  be  preferred. 
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Chapter  6 


Sensitivity  analyses  of  a 
metapopulation  model 

6.1  Introduction 

The  dynamics  of  most  populations  are  governed  by  interactions  between  demographic 
processes,  spatial  processes  (e.g.  dispersal),  and  stochastic  environmental  variability 
(e.g.  disturbance).  Each  of  these  processes  is  sufficiently  complex  that  ecologists  have 
tended,  depending  on  their  objectives,  to  focus  on  one  at  at  time.  For  instance,  con¬ 
servation  ecologists  tend  to  concentrate  on  demographic  structure  (e.g.  Mert.z,  1971; 
Crouse  et  ah,  1987;  Pascual  and  Adkison,  1994)  since  population  growth  or  decline  is 
of  primary  concern,  and  is  ultimately  determined  by  demographic  parameters.  Man¬ 
agers  of  harvested  populations,  such  as  fish  or  deer,  also  tend  to  focus  on  demography 
and  devise  management  strategies  that  are  based  on  parameters  such  as  such  as  sex, 
size,  or  age  (e.g.  Brenden  et  ah,  2007;  Collier  and  Krementz,  2007).  In  recent  years, 
interest  in  implementing  marine  protected  areas  has  shifted  some  attention  towards 
spatial  structure  (Roughgarden  and  Iwasa,  1986;  Thorrold  et  ah,  2001;  Fagan  and 
Lutscher,  2006).  Similarly,  growing  interest  in  metapopulation  theory  (Hanski,  1999) 
has  increased  the  number  of  studies  that  focus  on  spatial  processes  (Kritzer  and  Sale, 
2006). 

Ideally,  demographic  and  spatial  structure,  as  well  as  environmental  variability, 
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would  be  taken  into  account  when  deciding  how  best  to  manage  real  populations. 
In  addition,  since  management  efforts  are  inevitably  limited  by  resources  (i.e.  time, 
manpower,  money  etc.),  it  is  necessary  to  have  a  quantitative  way  of  choosing  the 
patches  and/or  stages  on  which  to  focus.  One  way  to  select  patches  is  according 
to  their  ability  to  contribute  to  metapopulation  growth  rate.  Of  the  many  different 
indices  used  for  population  assessment,  growth  rate  is  the  most  important  because  of 
the  power  of  exponential  growth  (Caswell,  2001,  Chap.  18). 

In  addition  to  demography  and  dispersal,  we  are  interested  in  how  patch  dynam¬ 
ics  are  influenced  by  the  presence  of  environmental  variability;  we  therefore  added 
stochastic  disturbance  to  our  model.  We  are  specifically  interested  in  how  elasticities 
to  individual  patch  demographic  parameters  change  with  relation  to  one  another.  We 
explore  the  case  where,  in  the  absence  of  disturbance,  one  patch  has  a  positive  growth 
rate  (the  “good”  patch)  and  the  other  has  a  negative  growth  rate  (the  “bad”  patch). 
We  use  the  “elasticity  ratio”  throughout  this  manuscript  to  quantitatively  assess  the 
relative  elasticities  of  growth  rates  in  the  good  and  bad  patches  to  their  respective 
parameters.  We  define  this  ratio  as  the  elasticity  of  growth  rate  to  parameters  in  the 
good  patch,  divided  by  the  elasticity  to  those  same  parameters  in  the  bad  patch. 

The  elasticity  ratio  (E)  has  the  potential  to  be  useful  when  deciding  on  the  best 
strategy  to  increase  metapopulation  growth  rate.  Assuming  all  things  are  equal  in 
the  good  and  bad  patch  (other  than  the  parameters  responsible  for  causing  one  to  be 
good  and  the  other  bad),  using  the  elasticity  ratio  is  straightforward.  In  this  simplest 
case,  if  E  >  1,  efforts  should  focus  on  the  good  patch,  and  if  E  <  1,  efforts  should 
focus  on  the  bad  patch.  Rarely  in  natural  metapopulations,  however,  are  all  things 
equal  among  patches.  For  instance,  it  may  be  more  costly  to  enact  changes  in  one 
patch  versus  another,  and  the  difference  in  cost  should  be  accounted  for  in  analyses. 
If  acting  in  patch  1  is  10  times  more  expensive  than  in  patch  2,  then  E  need  not  be 
less  than  1  but  only  less  than  0.10  to  suggest  efforts  focus  on  the  bad  patch.  Although 
we  examine  here  only  on  the  simplest  case  where  the  deciding  quantity  for  E  is  1,  it 
is  important  to  remember  that  the  decision-making  process  must  take  into  account 
the  relative  effect  of  managing  each  patch,  including  cost,  feasibility,  legality,  etc.. 
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In  this  manuscript,  we  combine  demographic  and  spatial  structure  in  a  matrix 
metapopulation  model.  We  are  not  the  first  to  explore  this  combination,  however. 
One  example  is  using  a  “megamatrix”  approach  (Pascarella  and  Horvitz,  1998;  Tul- 
japurkar  et  al.,  2003;  Petr,  2007),  introduced  by  Horvitz  and  Schemske  (1986).  In 
this  method,  both  the  spatial  and  demographic  dynamic  processes  are  included  in  the 
projection  matrix.  More  recently  Hunter  and  Caswell  (2005)  showed  how  to  combine 
demographic  and  spatial  structuring  in  a  way  that  simplifies  calculations,  particu¬ 
larly  calculations  involving  the  sensitivity  of  growth  rate  to  lower-level  parameters. 
Their  approach  uses  the  vec-permutation  matrix  to  rearrange  demographic  matrices 
arranged  by  patch  and  dispersal  matrices  arranged  by  stage,  thereby  allowing  both 
types  of  structure  within  the  same  model. 

Here,  we  build  on  the  work  of  Hunter  and  Caswell  (2005),  expanding  the  tech¬ 
nique  to  stochastic  models  and  calculating  the  elasticities  of  metapopulation  growth 
rate  to  individual  parameters  in  the  demographic  and  dispersal  matrices  (Caswell. 
2005,  2007).  Our  construction  facilitates  the  efficient  calculation  of  sensitivity  of 
metapopulation  growth  rate  to  demographic,  dispersal,  and  other  parameters.  The 
formulae  maintain  these  processes  as  distinct  matrices,  which  allows  for  more  intuitive 
exploration  of  the  dynamics  of  each  process. 

We  begin  by  examining  a  two-patch  metapopulation  with  no  disturbance  or  stage 
structure  using  a  deterministic  model.  We  then  determine  the  effects  of  stochastic 
disturbance  on  the  population  by  adding  a  variable  environment  to  the  model  and 
re-evaluating  our  results.  We  explore  the  elasticity  of  metapopulation  growth  rate 
to  migration  parameters  specifically  since  migration  dictates  the  demographic  con¬ 
nectivity  of  inetapopulation  patches  (Stacey  et  al.,  1997)  and  is  especially  critical  for 
long-term  inetapopulation  success  when  patches  are  subjected  to  stochastic  environ¬ 
mental  variability  (Howe  and  Davis,  1991;  Bascompte  et  al.,  2002;  Hill  et  al.,  2002). 
Finally,  we  include  stage  structure  to  more  closely  resemble  a  natural  population  and 
compare  the  outcomes  to  previous  simpler  versions  of  the  model.  Throughout  this 
study,  we  use  “patch  1”  to  refer  to  the  good  patch  and  “patch  2”  to  refer  to  the  bad 
patch. 
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6.2  One  Stage,  No  Disturbance 


Imagine  a  metapopulation  consisting  of  two  patches.  Let  rii(t)  be  the  population 
density  in  patch  i  at  time  t,  referred  to  hereafter  as  population  i.  Two  processes 
acting  sequentially  will  account  for  changes  in  First,  individuals  survive  and 

reproduce  with  a  net  per  capita  rate  Ri  in  patch  i.  Next,  a  proportion  mi  of  individuals 
from  patch  i  emigrate  to  patch  j.  Combining  these  two  processes,  we  can  write  the 
following  matrix  model  to  project  the  population  from  time  t  to  t  +  1: 

n(<  +  l)  =  An(t)  (6.1) 


where  n (t)  =  [ni(t) ,  n2(t)]J  and 


A  = 


Ri  (1  -  mi) 


Ri  mi 


R2m2 


\ 


R2  (1  -  m2) 


(6.2) 


In  our  model,  we  set  Ri  >  1  and  R2  <  1.  Thus  in  the  absence  of  migration,  population 
1  is  increasing  in  size  and  population  2  is  decreasing  in  size.  In  this  simple  model, 
the  metapopulation  growth  rate  A  is  the  dominant  eigenvalue  of  A,  and  the  stable 
patch  distribution  w  is  the  right  eigenvector  corresponding  to  A. 


We  used  sensitivity  analysis  to  determine  the  relative  effects  of  changes  to  R\  and 
R2  on  A.  The  sensitivities  of  A  to  the  elements  of  A  are  given  by  the  sensitivity 
matrix  with  entries 


s 


ij  ~ 


d\ 

O&ij 


(6.3) 


This  matrix  is 


S  = 


vw 


T 


v*w 


(6.4) 


where  v  is  the  left  eigenvector  corresponding  to  A  (Caswell,  2001).  The  entries  of  S 
are  the  sensitivity  of  A  to  the  corresponding  entry  in  A,  and  therefore  depend  on  a 
combination  of  the  parameters  that  make  up  that  entry.  A  more  interesting  analysis 
is  the  sensitivity  of  A  to  the  lower  level  parameters  that  compose  A.  From  Caswell 
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Figure  6.1:  Contour  plots  of  log  E :  red  contours  are  where  E  >  1  and  blue  contours 
are  where  E  <  1 .  Shaded  regions  are  where  A  >  1 . 


(2001,  Chap.  9),  we  calculate  the  proportional  sensitivity,  or  elasticity,  of  A  to  the 
lower  level  parameter  R,  as 


Ri  OX  Rj  sr  '  OX  0u,j 

X  ORi  X  Odii  OR , 

ij  j 

We  then  define  the  elasticity  ratio  of  A  to  Rl  as 

Ri  OX/OR , 

"  R2 ' 0X/0R2 


(6.5) 


(6.6) 


If  E  >  1,  then  A  is  more  elastic  to  changes  in  Ry.  If  E  <  1,  then  A  is  more  elastic 
to  changes  in  R2.  The  value  of  E  can  be  used  to  choose  patches  on  which  to  focus 
efforts  when  the  goal  is  to  increase  metapopulation  growth  rate  efficiently. 
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Figure  G.l  shows  log  E  for  different  combinations  of  net  reproductive  rates,  with 
rni  and  ra2  ranging  from  zero  to  one.  When  migration  out  of  the  good  patch  was  low, 
log  E  tended  to  be  less  than  0.  (Fig.  6.1).  That  is,  if  individuals  in  the  good  patch 
tend  to  stay  in  that  patch,  then  regardless  of  the  fraction  of  individuals  who  migrate 
out  of  the  bad  patch  (m2),  metapopulation  growth  rate  is  more  sensitive  to  changes 
in  Ri  than  /?2.  As  migration  out  of  the  good  patch  mi  increases,  the  proportion  of 
individuals  located  in  the  bad  patch  also  increases  and  E  becomes  less  than  1. 

Using  equations  (6.4)-(6.6),  we  can  write  an  expression  for  E  in  terms  of  the  model 
parameters  Ri  and  m*  and  entries  of  the  right  and  left  eigenvectors  v  and  w: 

E  =  RjWi  K  (1  ~  mQ  +  u2mi]  ^  ^ 

R2  w-2  [v2  (1  -  m2)  +  mm2] ' 

The  equation  for  the  line  where  E  =  1  (the  dashed  line  in  Fig.  6.1)  is 

m2  =  1  +  (mL  -  1)  •  (6-8) 

n-2 

To  the  left  of  this  line  (as  mi  decreases),  E  >  1  and  increasing  R\  results  in  the 
largest  proportional  increase  in  A.  To  the  right  of  this  line  (as  mi  increases),  E  <  1 
and  increasing  /?2  results  in  the  largest  proportional  increase  in  A. 

The  values  for  A  are  also  of  importance  for  patch  prioritization,  as  the  transition 
from  A  <  1  to  A  >  1  indicates  where  the  metapopulation  switches  from  declining 
to  growing.  As  one  might  expect,  increasing  net  reproductive  rate  in  either  patch 
increases  the  area  of  parameter  space  where  A  >  1  (shaded  areas  in  Fig.  6.1). 

An  unexpected  result  is  that  there  are  instances  when  E  <  1  and  A  >  1  (Fig.  6.1B, 
6. ID).  That  is,  there  are  times  when  the  metapopulation  growth  rate  is  positive  and  is 
more  sensitive  to  net  reproductive  rate  in  the  bad  patch.  In  these  instances,  migration 
rates  are  such  that  the  majority  of  the  metapopulation  is  found  in  the  bad  patch. 
This  can  be  seen  if  one  imagines  a  1:1  line  in  any  of  the  plots  of  Fig.  6.1;  this  line 
is  where  mi  =  m2.  Below  the  line  there  are  more  individuals  in  patch  2,  and  above 
this  line,  there  are  more  individuals  in  patch  1.  Cases  where  E  <  1  and  A  >  1  are 
always  below,  suggesting  that  a  necessary  condition  for  this  result  is  that  the  majority 
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Patch  2  )  R,(l-m,)<\(t) 


R2»i2  6,(t) 


Figure  6.2:  Life  cycle  graph  for  one-stage,  two-patch  model  with  stochastic  distur¬ 
bance. 

of  the  metapopulation  is  found  in  the  bad  patch.  One  might  assume  that  if  the 
majority  of  individuals  are  in  the  bad  patch  then  growth  rate  of  the  metapopulation 
would  be  negative;  however,  this  is  clearly  not  the  case.  This  result  indicates  that 
increasing  /?2  has  a  larger  effect  because  it  impacts  more  individuals  than  increasing 
R[.  It  also  suggests  that  not  only  should  individual  patch  growth  rates  be  considered 
when  determining  a  patch’s  “goodness”  or  "badness”,  but  also  the  distribution  of 
individuals  among  patches  within  the  mctapopulation. 

6.3  One  Stage,  With  Disturbance 

Metapopulations  in  nature  are  subject  to  stochastic  environmental  variability.  We 
would  therefore  like  to  add  stochastic  disturbance  to  our  model  and  reassess  the 
elasticity  ratio.  One  typical  way  that  disturbance  affects  a  metapopulation  is  by 
reducing  its  size.  Here  we  model  disturbance  as  a  random  event  with  probability  p 
that  acts  by  reducing  population  size  in  patch  i  by  a  proportion  6t(t)  (Fig.  6.2). 

To  construct  6(t)  =  [<Ji(£),<52(£)]t.  we  first  let  be  an  indicator  variable  for  the 
event  that  population  i  is  disturbed,  such  that 


Xi(t)  = 


1  if  population  i  is  disturbed 


(6.9) 


0  otherwise 
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We  then  collect  these  events  into  the  random  vector  x(t),  and  set 


<5(0  =  1  -  Dx(t),  D  €  [0,1]. 


(6.10) 


Thus  D  is  a  measure  of  disturbance  intensity;  given  that,  patch  i  is  disturbed,  then 
Xi(t)  =  1  and  the  population  in  that  patch  is  reduced  by  the  proportion  D.  S(t)  is 
therefore  the  number  that  survive  the  disturbance,  equal  to  1  —  D.  Large  values  of 
D  result  in  high  disturbance  intensity. 

We  next  assume  that  x(t)  is  drawn  from  a  bivariate  Bernoulli  distribution  (Mar¬ 
shall  and  Olkin,  1985).  We  set  the  expectation 


E(x(t)|  = 


(6.11) 


so  that  the  patches  are  disturbed  with  equal  probability  at  each  time.  One  can 
imagine  that  the  occurrence  of  disturbance  at  two  patches  may  be  independent,  or 
it  may  be  positively  or  negatively  correlated.  If.  for  instance,  the  disturbance  event 
is  a  large  weather  event  such  as  a  hurricane,  it  might  affect,  all  of  the  patches  in  a 
metapopulation  and  therefore  covariance  of  patch  disturbances  (c)  would  be  positive. 
With  c  =  cov(xi,X2),  the  variance-covariance  matrix  for  x(t)  is  the  constant  matrix 


var[x(f)]  = 


p(l-p) 


P(1  ~P) 


(6.12) 


One  can  show  (Appendix)  that  p  and  c  must  satisfy  the  inequalities 


< 

p{  1  -  p ), 

(6.13) 

> 

-p\ 

(6.14) 

> 

-(1  -p)‘\ 

(6.15) 

and  that  the  sum  of  probabilities  of  all  possible  disturbance  events  equals  1.  Inequal¬ 
ities  (6.13)-(6.15)  define  a  two-dimensional  parameter  space  of  all  allowable  combina- 
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Figure  6.3:  Allowable  combinations  (shaded),  as  defined  by  inequalities  (6.13)-(6.15). 
of  the  probability  of  disturbance  {p)  and  the  covariance  of  disturbance  events  (c)  at 
the  two  patches  in  model  (6.9)-(6.16). 


tions  of  p  and  c  (Fig.  6.3). 


Including  disturbance,  the  deterministic  projection  matrix  (6.2)  becomes  the  stochas¬ 
tic  matrix 

R]  ( 1  -  mi )  S\  ( t )  R‘2  m2  Si  ( t ) 

Ri  rrti  S2(t)  R2{  1  -m2)S2(t) 


A*  — 

The  stochastic  growth  rate 


\ 

/ 


(6.16) 


log  A.,  =  Jim  —  log  ||  A/-!  •••  A0n,)|| 
/  — >00  1 


(6.17) 


is  (except  under  bizarre  circumstances)  the  long-term  average  growth  rate  of  every 
realization  of  the  model  with  probability  1  (Furstenberg  and  Kesten.  1960;  Cohen. 
1976:  Tuljapurkar  and  Orzack.  1980:  Caswell,  2001).  The  elasticity  of  the  stochastic 
growth  rate  to  the  net  reproductive  rate  in  patch  i  is 


d  log  A,  J_  ^4  ftjVT(<  +  l)j^w(Q 

d log  Ri  T^x  T  r(t)\T(t  +  l)w(f  +  1) 


(6.18) 


(Caswell.  2005.  pg.  80).  Here  v(f)  and  w (t)  are  the  stochastic  analogs  to  the  left 
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and  right  eigenvectors  of  the  deterministic  projection  matrix,  and  r(t)  is  the  one-step 
growth  rate 


_  l|Atw(Q|| 

llw(*)ll 


(6.19) 


Equation  (6.18)  gives  the  proportional  change  in  log  A,,  with  net  reproductive  rate 
over  many  stochastic  simulations.  Bigger  values  for  (6.18)  indicate  bigger  changes  in 
log  A.,  with  a  proportional  change  in  /?,.  We  define  the  stochastic  elasticity  ratio  of 


log  As  to  the  Ri  as 


E,  = 


<9  log  As/51og  R  i 


(6.20) 


<9  log  As/<91og  R2 

As  in  the  deterministic  case,  if  Es  >  1  then  log  A.,  is  more  elastic  to  changes  in  R\ . 
Conversely,  if  E,  <  1,  then  log  A*  is  more  elastic  to  changes  in  R2. 


We  calculated  log  As  and  the  elasticity  ratio  (6.20)  for  R i  —  2.5,  R2  —  0.9,  D  —  0.9, 
and  for  four  migration  scenarios  in  which  migration  is  high  or  low'  from  the  patches 
and  either  equal  or  not: 


(mi, m2)  €  < 


(0.1,  0.1) 
(0.1,  0.9) 
(0.9,  0.1) 
(0.9,  0.9) 


(6.21) 


As  in  the  deterministic  case,  the  elasticity  ratio  is  correlated  to  the  distribution  of 
individuals  among  patches  in  the  metapopulation.  Determining  where  the  majority  of 
individuals  are  located  is  not  as  straightforward  in  the  stochastic:  case.  Rather  than 
determine  this  based  on  migration  rates,  we  must  instead  calculate  the  long-term 
average  patch  distribution. 


In  (6. 4 A),  there  are,  on  average,  more  individuals  in  patch  1  than  in  patch  2  since 
few  individuals  leave  patch  1  mi  =  0.1  and  most  individuals  from  patch  2  migrate 
to  patch  1  m2  =  0.9.  As  expected,  elasticity  ratio  is  greater  than  1  and  changes 
to  the  good  patch  will  result  in  proportionally  greater  increases  in  metapopulation 
growth  rate.  In  the  case  where  migration  rates  are  equal  (6.4B-C),  individuals  are, 
on  average,  distributed  equally  between  the  two  patches.  Again,  the  elasticity  ratio 
is  greater  than  1. 
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Figure  6.4:  Contour  plots  of  log  Es  over  the  range  of  possible  values  for  disturbance 
parameters  p  and  c,  for  four  combinations  of  mi  and  ra2.  Shaded  regions  indicate 
where  log  A*  >  0.  Blue  contours  are  where  Es  <  1:  red  contours  are  where  Es  >  1. 
For  all  panels,  R\  2.5,  R2  =  0.9,  and  D  =  0.9.  Here,  contours  are  log  Es  to  show 
the  details  of  the  response  of  log  A.s  at  small  values  and  for  comparing  results  to  the 
no  disturbance  case  (Fig.  6.1). 
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The  result  changes  if  emigration  from  patch  1  is  high  and  emigration  from  patch  2 
is  low  (Fig.  6.4D).  Migration  rates  in  this  case  result  in  more  individuals,  on  average, 
in  patch  2.  This  is  because  migration  out  of  patch  1  is  high  (rri\  =  0.9  and  most 
individuals  from  patch  2  do  not  migrate  (m^  =  0.1).  Consequently  Es  <  1  for 
all  combinations  of  disturbance  parameters  p  and  c.  These  results  parallel  those  of 
the  deterministic  model:  Es  <  1  only  if  there  are  more  individuals  on  average  in 
population  2  than  in  population  1.  Changes  to  parameters  in  patch  2  will  affect 
more  individuals  of  the  metapopulation,  resulting  in  proportionally  greater  increase 
in  metapopulation  growth  rate. 

For  the  reproductive  rate  and  disturbance  intensity  used  for  Fig.  6.4,  the  stochastic 
metapopulation  growth  rate  is  positive  only  when  the  probability  of  disturbance  is  low 
(shaded  areas,  Fig.  6.4).  In  addition,  the  probability  of  disturbance  appears  to  have 
a  greater  effect  on  population  growth  rate  than  the  covariance  of  patch  disturbance. 
The  set  of  disturbance  parameters  p  and  c  for  which  log  As  >  0  is  largest  when  on 
average  population  1  is  larger  than  population  2.  The  set  shrinks  as  individuals 
become  more  evenly  distributed  between  the  two  patches.  When  population  2  is 
larger,  as  in  Fig.  6.4D,  log  A.,  >  0  only  for  the  smallest,  values  of  the  disturbance 
probability  p. 

Although  there  is  only  a  small  amount  of  variability  in  Es  within  a  given  panel. 
Ea  varies  dramatically  among  the  four  panels  (note  the  log  scale).  Migration  rates 
appear  to  affect  Es  to  a  much  greater  extent  than  disturbance  parameters.  Es  <  1 
only  when  emigration  was  high  from  patch  1  and  low  from  patch  2  (Fig.  6.4D),  and 
in  that  migration  scenario  log  A.,  <  0  for  most  combinations  of  p  and  c.  There  is  a 
small  set  of  disturbance  parameter  values  (p  <  0.05,  c  ~  0)  where  both  log  As  >  0 
and  Es  <  1;  this  occurs  at  the  lowest  values  of  p.  when  disturbance  is  so  improbable 
that  results  from  the  stochastic  model  are  comparable  to  those  of  the  deterministic 
model  (Fig.  6.1B). 

To  further  understand  the  effects  of  migration  on  metapopulation  dynamics,  we 
calculated  the  elasticity  of  log  As  to  both  mi  and  m-2  under  the  same  four  migration 
scenarios  explored  in  Fig.  6.4  (Fig.  6.5).  We  did  not  calculate  elasticity  ratios  since 


114 


Elasticity  of  log  A.  to  ml  Elasticity  of  log  Ato  /», 
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m2  =  0.9 
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Figure  6.5:  Elasticity  of  log  A.s  to  changes  in  migration  for  four  migration  scenarios. 
Blue  contours  are  where  elasticity  is  negative:  red  contours  are  where  elasticity  is 
positive.  For  all  panels.  H\  =  2.5,  R>  —  0.9,  and  D  =  0.9. 


115 


elasticity  values  were  both  positive  and  negative,  thereby  complicating  the  interpreta¬ 
tion  of  such  a  ratio.  In  general,  when  migration  was  such  that  there  were,  on  average, 
an  equal  number  of  individuals  in  the  two  patches  (6.5C-6.5F)  or  there  were  more 
individuals  in  patch  2  (6.5G-6.5H),  elasticity  to  was  less  than  0  and  elasticity  to 
m-2  was  greater  than  zero.  These  results  indicate  that  log  As  was  negatively  impacted 
by  increasing  migration  from  the  good  patch,  and  positively  impacted  by  increasing 
migration  from  the  bad  patch.  Results  were  more  complicated  when  migration  re¬ 
sulted  in  more  individuals  on  average  in  patch  1  (6.5A-6.5B).  As  covariance  decreased, 
elasticity  to  mi  switched  from  negative  to  positive  and  elasticity  to  m-2  switched  from 
positive  to  negative.  This  corresponds  to  the  good  and  bad  patches  essentially  switch¬ 
ing  their  "goodness”  and  “badness”  when  covariance  is  negative,  so  that  increasing 
migration  from  patch  1  results  in  increased  metapopulation  growth  rate,  and  increas¬ 
ing  migration  from  patch  2  causes  decreases  in  metapopulation  growth  rate. 


6.4  Two  Stages,  With  Disturbance 

The  softshell  clam,  Mya  arenaria ,  is  a  commercially  important  bivalve  commonly 
found  in  New  England  estuaries.  M.  arenaria' s  life  cycle  is  typical  of  nearshore  marine 
benthic  invertebrates  (Thorson,  1950;  Abraham  and  Dillon,  1986).  It  is  characterized 
by  a  relatively  sedentary  adult  stage,  with  adults  highly  aggregated  into  patches  of 
suitable  habitat.  Dispersal  between  these  populations  is  accomplished  by  a  short¬ 
lived  larval  stage.  In  many  species,  larvae  are  produced  in  vast  quantities  during  a 
short  reproductive  season.  Almost  all  of  the  larvae  die  before  recruiting  to  the  adult 
phase. 

To  model  such  a  population  requires  at  least  two  patches,  with  dispersal  between 
them,  and  two  stages:  adults  and  larvae/juveniles  (Fig.  6.6).  We  assume  that  demog¬ 
raphy.  migration,  and  disturbance  act  sequentially  (in  the  intervals  (Mi),  (tuh)  and 
(M  t  +  1)  respectively)  within  a  given  projection  interval  (M  +  1). 

Between  times  t  and  ti,  the  population  undergoes  demographic  processes.  The 
demographic  transitions  within  population  i  are  given  by  the  matrix  B,,  with  entries 
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t  + 


Figure  6.6:  Life  cycle  graph  for  two-stage,  two-patch  model  with  stochastic  distur¬ 
bance.  B.  M.  and  D  denote  the  block-diagonal  matrices  for  demography,  migration, 
and  disturbance,  respectively  (Equations  (6.25),  (6.29),  and  (6.31)). 
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composed  of  adult  survival  (cr*),  per  capita  larval  production  and  survival  ($ ),  and 
survival  and  maturation  of  new  recruits  (7*): 


Bj  = 


(6.22) 


If  we  define  the  entries  of  the  vector  n(t)  to  be  the  number  of  individuals  in 

stage  j  of  population  i  at  time  t  and  arrange  the  elements  as 


n  (t)  = 


*  nn(t)  ^ 
nn(t) 
n2i{t) 


^22  (t)  J 


then  the  demographic  transitions  are  described  by 


(6.23) 


n(tj)  =  Bn(t) 


(6.24) 


where 


(6.25) 


Next,  individuals  migrate  between  patches.  Let  Mj  be  the  matrix  of  migration 
rates  for  individuals  in  stage  j.  A  simple  model  for  migration  has  these  stage- j 
individuals  migrating  from  patch  i  at  the  per  capita  rate  rriij.  Thus, 


Mj 


1  —  rriij  m<2j 

mij  1  —  m2j 


(6.26) 
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If  we  now  rearrange  the  vector  n (t\)  as 


then 


n(fj) 


'  nn(^i)  ^ 
n2]{ti) 
nn{t\ ) 
y  ^22(^1)  j 


n(t2)  =  Mn(ti), 


(6.27) 


(6.28) 


where 


(6.29) 


Since  M.  arenaiia  adults  are  sedentary,  we  set  M2  =  I.  To  simplify  notation,  we 
will  set  mn  =  m  1  and  m2 1  =  m2  from  here  on. 

Finally,  disturbance  reduces  the  number  of  individuals  in  population  i  by  the 
fraction  Si,  as  described  in  equations  (6.10)  through  (6.15).  If  we  return  to  the 
original  arrangement  (6.23),  we  then  have 


n  (t  +1)  —  KVn(f2). 


(6.30) 


where 


and 


Bt  = 


D  i(t)  0 
0  D2(*) 


D  i(t) 


8i(t)  0 

0  Si(t) 


(6.31) 


(6.32) 


To  convert  between  the  vectors  n  and  n.  conversions  that  are  required  at  every 
time  step,  we  employ  the  vec-permutation  matrix  P  (Henderson  and  Searle,  1981): 


n  =  Pn, 
=  PTn. 


n 
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(6.33) 

(6.34) 


For  the  two-patch,  two-stage  case 


P  = 


1 

0 

0 

0 


0  0  0  ^ 
0  1  0 
1  0  0 
0  0  1  ) 


(6.35) 


At  last,  combining  the  demographic,  migration,  and  disturbance  processes  gives 
the  model 

n(t  +  1)  =  Atn(i)  (6.36) 


where 

At  =  DtPTMP®. 


(6.37) 


This  decomposition  of  the  projection  matrix  into  its  demographic,  migratory  and 
disturbance  components  is  useful.  It  allows  us  to  use  matrix  calculus  to  calculate  the 
elasticity  of  the  stochastic  metapopulation  growth  rate  to  each  of  the  demographic 
parameters  (Caswell,  2007). 

The  elasticity  of  log  As  to  parameter  changes  is 


<91ogAs  1  >— -i  0(f)  wT(f)  ®  vT(f  +  1)  <9  vecA< 

<91og0T  T-,°°  T  ^  r(t)  \T{t  +  1)  w(t  +  1)  d6J 


where  0  is  a  column  vector  of  parameters,  and  the  operator  vec(-)  stacks  the  columns 
of  a  matrix  on  top  of  each  other  (Caswell,  unpublished).  The  matrix  dvecAt/d6r  is 
the  derivative  of  the  projection  matrix  at  each  time  step  with  respect  to  each  lower- 
level  parameter  in  6.  In  our  model,  At  is  the  product  of  several  matrices,  therefore 
the  calculation  of  dvecAt/dO1  is  not  trivial. 

Taking  the  differential  of  both  sides  of  equation  (6.37)  gives 


dA  =  (dD)  PtMPB  +  DPT(dM)  PB  +  OPTMP  (dM) .  (6.39) 


Multiplying  the  first  and  third  terms  in  equation  (6.39)  by  the  identity  matrix  leaves 


120 


those  terms  unchanged. 


dA  =  I  (dD)  PtMP1  +  DPT(dM)  P®  +  DPTMP  (dB)  I.  (6.40) 

We  can  then  apply  the  vec  operator  to  each  term  in  (6.40)  .  Using  the  fact  that 
vec(ABC)  =  (CT  0  A)  vecB,  we  obtain  an  equation  for  dve cA  in  terms  of  the  com¬ 
ponent  matrices: 


dve  cA 


[(PtMPB)  ®  I]  dvecD 


+ 


(PB)t0  (DPt) 


dvecM 


+  [IT  0  (DPtMP)]  dvecB. 


(6.41) 


Using  the  chain  rule,  along  with  the  first  identification  theorem  of  Magnus  and 
Neudecker  (1985)  then  gives 


dve  cA 


[(PtMPB)  01] 


+ 


dvecD 

d0T 

(PB)T  0  (dpt) 


1  dvecM 
dOr 


+  [IT  0  (dptmp)] 


dvecB 

d0J 


(6.42) 


The  matrices  dvecD/d^1 .  dvecM/d0T  and  dvecB/dt?1  can  be  rewritten  as  in  terms 
of  their  component  matrices  (i.e.  D,.  M,.  and  B,-).  For  example. 


dvecD 

d.0T 


2 

X)(H®I2) 

2=0 


dvecD  i 


(6.43) 


where 


H  =  (I2  0  P)  (vecEj  0  I2)  (6.44) 

and  E,  is  the  2x2  matrix  with  every  entry  zero,  save  the  ( i ,  i)-th,  which  is  1  (Magnus 
and  Neudecker,  1985).  To  calculate  dvecM/d0T,  simply  replace  D  and  D,  in  (6.43) 
by  M  and  M,  respectively.  Analogous  steps  yield  dvecB/d#7.1 

'The  procedure  described  in  equations  (6.22)  through  (6.44)  are  valid  for  the  two-stage,  two- 
patch  case.  The  elasticity  equation  corresponding  to  (6.42)  for  the  general  s-stage,  A-pat.ch  case  is 
more  complicated  but  can  be  similarly  derived. 
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We  analyzed  model  (6.36)-(6.37)  with  the  parameters  set  at:  cq  =  0.9,  a 2  =  0.3, 
71  =  0.8,  72  =  0.24,  =  5.6,  and  /?2  =  7.5.  We  chose  these  parameter  values  so 

that  (1)  the  individual  patch  growth  rates  would  be  similar  to  the  net  reproductive 
rates  for  the  one-stage  case  above,  i.e.  R\  =  2.5  and  i?2  =  0.9;  and  (2),  so  that 
they  would  roughly  comport  with  the  results  of  Ripley  and  Caswell  (2006)  who  esti¬ 
mated  demographic  parameters  for  M.  arenaria  from  field  studies.  Migration  rates 
are  difficult  to  obtain  for  M.  arenaria ,  as  they  are  for  all  benthic  invertebrates  with 
pelagic  larvae,  but  the  values  m\  =  m2  =  0.1  are  arguably  realistic.  We  assumed 
that  disturbance  affects  all  patches  and  stages  with  intensity  D  =  0.9,  and  set  the 
probability  of  disturbance  at  the  two  patches  to  p  =  0.5. 

When  the  covariance  of  patch  disturbance  c  is  positive,  both  patches  tend  to  be 
disturbed  simultaneously.  This  preserves  the  inherently  higher  quality  of  population 
1,  and  the  resulting  elasticities  of  logAs  to  population  1  parameters  (Fig.  6.7A)  are 
much  larger  than  the  elasticities  to  population  2  parameters  (i.e.  Es  >>  1).  When 
covariance  is  negative,  the  patches  tend  to  be  disturbed  at  different  times.  As  a  result, 
population  2  can  be  temporarily  of  higher  quality  than  population  1.  The  elasticities 
of  log  As  to  population  2  parameters  are  therefore  larger  than  they  are  when  c  >  0 
(Fig.  6.7B).  Nevertheless  elasticities  to  population  2  parameters  are  still  much  smaller 
than  those  to  population  1  parameters.  These  results  parallel  those  we  obtained  for 
the  one-stage  model  (cf.  Fig.  6.4C).  In  both  disturbance  scenarios,  the  elasticity  of 
logAs  was  negative  to  m\  and  positive  to  m2;  the  magnitude  of  these  elasticities  was 
quite  small,  however,  compared  to  the  other  parameters  investigated. 

We  also  explored  the  effects  of  dispersal  scenarios  that  resulted  in  the  majority 
of  the  metapopulation  density  being  in  patch  1  ( mi  =  0.1,  m2  =  0.9;  Fig.  6.8A) 
or  patch  2  (mi  =  0.9,  m2  =  0.1;  Fig.  6.8B).  We  chose  a  low  probability  of  intense 
disturbance  ( p  =  0.15,  D  =  0.9)  and  no  covariance  between  patches  (c  =  0).  In  both 
dispersal  scenarios,  log  As  is  more  elastic  to  changes  in  population  1  parameters  than 
it  is  to  the  corresponding  parameters  in  population  2.  When  the  population  in  patch 
1  is  larger  than  the  population  in  patch  2  (Fig.  6.8A),  elasticities  to  population  2 
parameters  are  much  smaller  than  to  population  1.  When  the  majority  of  individu- 
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Figure  6.7:  Elasticity  of  log  Xs  to  changes  in  lower-level  parameters.  0  under  two  dis¬ 
turbance  scenarios:  probability  of  disturbance  p  is  0.5,  with  covariance  either  positive 
(A)  or  negative  (B).  White  bars  are  elasticities  in  patch  1;  gray  bars  are  elasticities 
in  patch  2.  In  both  plots,  mi  =  0.1,  m2  —  0.1,  D  =  0.9,  and  logAs  <  0  (-0.27  and 
-0.23  for  (A)  and  (B)  respectively). 


als  are  found  in  population  2  (Fig.  6.8B),  elasticity  to  population  2  parameters  are 
much  closer  to  those  to  population  1,  but  logAs  still  is  most  elastic  to  population  1 
parameters.  For  both  migration  scenarios,  the  elasticity  of  log  A*  was  negative  to  m\ 
and  positive  to  m2.  The  magnitudes  of  these  elasticities  were  quite  small  for  t  he  mi¬ 
gration  scenario  that  resulted  in  more  individuals  in  patch  1  (6.8A),  and  much  larger 
for  the  second  migration  scenario  that  resulted  in  more  individuals  in  patch  2  (6.8B). 
When  the  model  complexity  was  increased  by  adding  a  second  stage,  the  elasticities 
of  stochastic  metapopulat  ion  growth  rate  to  population  1  parameters  increased,  while 
those  same  elasticities  to  population  2  parameters  deereased.  Even  when  the  major¬ 
ity  of  individuals  were  in  the  bad  patch,  growth  rate  elasticity  was  still  greatest  to 
parameters  in  the  good  patch. 


6.5  Discussion  &  Conclusions 


We  have  demonstrated  how  to  incorporate  stochastic  disturbance  into  matrix  popu¬ 
lation  model  with  two  patches.  Our  results  suggest  that  stochastic  metapopulation 
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Figure  6.8:  Elasticity  of  logAs  to  changes  in  lower-level  parameters  (0)  under  two 
migration  scenarios:  (A)rai  small  and  m2  large,  resulting  in  more  individuals  in 
patch  1,  and  (B)mi  large  and  m2  small,  resulting  in  more  individuals  in  patch  2. 
White  bars  are  elasticities  in  patch  1;  gray  bars  are  elasticities  in  patch  2.  In  both 
plots,  p  =  0.15,  c  =  0,  and  D  =  0.9,  and  logAs  >  0  (0.53  and  0.05  for  (A)  and  (B) 
respectively). 


growth  rate  tends  to  be  more  sensitive  to  the  population  that  holds  the  larger  propor¬ 
tion  of  the  total  metapopulation.  This  is  true  with  or  without  stochastic  disturbance, 
although  as  expected  metapopulation  growth  rate  is  lowrer  when  disturbance  occurs. 
When  a  second  stage  was  added,  the  elasticities  of  stochastic  metapopulation  growth 
rate  to  population  1  parameters  increased,  while  those  same  elasticities  to  population 
2  parameters  decreased.  Metapopulation  growth  rate  is  most  elastic  to  adult  survival 
in  population  1  for  all  scenarios  we  examined.  If  the  majority  of  the  metapopula¬ 
tion  is  located  in  the  bad  patch,  then  the  elasticity  to  changes  in  parameters  of  that 
population  increase  but  do  not  surpass  elasticity  to  changes  in  the  good  patch. 

Previous  studies  have  reported  that  the  distribution  of  individuals  is  important 
for  determining  the  relative  effects  of  good  and  bad  patches,  but  have  approached 
the  problem  from  a  different  perspective.  Managers  sometimes  use  abundance  as 
an  indicator  of  habitat  quality  (i.e.  for  determining  good  and  bad  patches)  because 
it  is  easily  determined  by  surveys  (e.g.  Peres,  2001).  This  strategy  is  not  valid, 
however,  since  low-abundance  good  patches  can  support  high-abundance  bad  patches 
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via  dispersal  (Pulliam,  1988).  When  managers  neglect  to  identify  processes  that 
drive  changes  in  abundance,  the  may  incorrectly  infer  the  status  of  patches  based 
solely  on  abundance  (Van  Horne,  1983).  Our  results  suggest  that  the  problem  can  be 
reversed;  one  may  assume  it  ineffective  to  invest  in  improving  the  bad  patch,  while  if 
abundances  are  factored  in,  focusing  management  efforts  on  the  bad  patch  will  result 
in  the  most  efficient  increases  in  metapopulation  growth  rate. 

A  more  common  way  of  referring  to  good  and  bad  patches  is  as  “sources”  and 
“sinks”.  The  source/sink  literature  (Runge  et.  ah,  2006;  Howe  and  Davis,  1991:  Pul¬ 
liam,  1988)  and  intuition  suggest  that  if  one  must  choose  between  focusing  manage¬ 
ment  efforts  on  a  source  or  a  sink,  one  should  always  choose  the  source.  Our  results 
comport  with  this  choice,  especially  when  there  is  no  available  information  about 
exchange  rates  via  migration.  This  result  does  not,  however,  preclude  the  necessity 
of  the  sink  populations  to  long-term  metapopulation  persistence.  In  a  constant  en¬ 
vironment  with  density-dependent  recruitment,  individuals  unable  to  settle  at  source 
populations  due  to  high  densities  can  migrate  to  sinks.  In  this  wav,  the  presence 
of  a  sink  results  in  larger  overall  metapopulation  size  (Pulliam,  1988).  Emigration 
from  the  source  patch  becomes  beneficial,  as  it  offers  “insurance”  against  catastrophe 
(Levin  et  al.,  1984).  This  is  especially  true  if  the  environmental  variability  is  spatially 
negatively  correlated  (Wiener  and  Tuljapurkar,  1994).  We  found  that  this  was  true 
for  both  one-stage  and  two-stage  cases  we  examined.  As  covariance  of  patch  distur¬ 
bances  decreases,  the  elasticity  ratio  also  decreases.  That  is,  effects  of  parameters 
in  the  bad  patch  increase,  and  effects  of  those  same  parameters  in  the  good  patch 
decrease  (Figs.  6.4,  6. 7-6. 8). 

The  benefits  of  sinks  are  even  more  pronounced  when  the  environment  varies 
over  time.  If  population  density  in  the  source  falls  below  a  sustainable  level  due 
to  a  stochastic  event,  the  sink  can  serve  as  a  refuge  for  the  metapopulation  and 
provide  emigrants  to  recolonize  the  source  (Runge  et  al.,  2006).  Based  on  their  model 
of  source/sink  dynamics,  Howe  and  Davis  (1991)  similarly  concluded  that  although 
sinks  may  not  persist  independently,  they  contribute  to  metapopulation  size  and 
longevity.  Runge  et  al.  (2006)  also  found  that  sinks  may  be  critical  to  metapopulation 
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persistence.  They  found  that  although  populations  may  appear  to  be  a  sink  based  on 
demographic  parameters,  there  may  be  a  sufficient  number  of  emigrants  to  classify 
it  as  a  source.  We  similarly  found  that  migration  rates  to  and  from  the  bad  patch 
are  critical  for  determining  their  roles  in  affecting  stochastic  metapopulation  growth 
rate.  Migration  rates  that  result  in  more  individuals  in  the  bad  patch  also  result  in 
decreases  in  the  elasticity  ratio. 

Migration  strongly  affects  metapopulation  growth  rate,  in  some  cases  surpassing 
in  magnitude  the  elasticity  to  demographic  parameters  (e.g.  Fig.  6.8B).  The  direc¬ 
tional  relationship  between  elasticity  and  metapopulation  growth  rate  is  dictated  by 
patch  quality.  Elasticity  of  metapopulation  growth  rate  to  migration  from  the  high- 
quality  patch  tended  to  be  negative,  while  the  same  elasticity  for  the  low-quality  patch 
was  positive.  These  relationships  switched  when  covariance  of  patch  disturbance  was 
negative,  further  indicating  that  the  interactive  effects  of  dispersal  and  demogra¬ 
phy  cannot  be  ignored.  Results  here  strongly  suggest  that  migration  rates  play  a 
substantial  role  in  patch  dynamics  and  metapopulation  growth  rate,  and  therefore 
offer  further  evidence  that,  metapopulation  studies  should  move  towards  obtaining 
estimates  of  connectivity  to  advance  the  field. 

Although  the  elasticity  ratio  is  a  useful  metric  for  determining  the  relative  ef¬ 
fects  of  patches  on  metapopulation  growth  rate,  other  factors  might,  be  necessary  to 
consider  when  deciding  where  to  focus  management  efforts.  If  one  patch  is  inside 
a  protected  area  (e.g.  a  national  park),  it  may  be  illegal  to  take  steps  towards  in¬ 
creasing  parameters  such  as  survival  or  reproduction.  One  patch  may  be  much  more 
accessible  than  another  patch  due  to  geography,  making  it  more  realistic  for  managers 
to  focus  on  that  patch  irrespective  of  the  elasticity  ratio.  Other  costs  associated  with 
management,  efforts  might  result  in  similar  situations  where,  once  “real  world”  fac¬ 
tors  are  considered,  the  outcome  of  the  elasticity  ratio  may  be  a  poor  indicator  of 
the  most  logical  management  plan.  This  ratio  is  therefore  useful  as  a  tool  for  guiding 
management  decisions,  along  with  associated  costs,  legality  issues,  accessibility,  and 
other  considerations. 
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6.6  Appendix 


6.6.1  Bivariate  Bernoulli  Distribution 

There  are  four  possible  disturbance  scenarios: 

1.  Both  populations  are  disturbed  (6^62  =  d ); 

2.  Population  1  is  disturbed  and  Population  2  is  not  (5]  =  d,S2  =  1); 

3.  Population  2  is  disturbed  and  Population  1  is  not  (5j  =  1,  62  =  d); 

4.  Both  populations  are  not  disturbed  (<5i,52  =  1). 

If  we  assign  each  of  these  four  possibilities  a  probability  Qj,  we  can  establish  a  set  of 
rules  under  which  disturbance  acts.  First, 

9i  +  92  +  93  +  94  =  1  (6.45) 

where  0  <  qj  <  1.  That  is,  the  four  disturbance  scenarios  described  above  are  the  only 
four  possible  scenarios  that  can  occur.  Second,  disturbance  must  be  equally  probable 
in  the  two  populations  so  that  any  differences  in  sensitivity  of  metapopulation  log  A*  to 
the  demography  of  the  populations  is  due  to  their  inherent  “goodness”  and  “badness”, 
i.e.  their  rather  than  the  nature  of  the  disturbance.  If  this  is  true,  then 


9i  +  92  =  P 

(6.46) 

9i  +  93  =  P 

(6.47) 

and  therefore 

(6.48) 

92  =  93 

(6.49) 

Finally,  we  establish  the  covariance  of  disturbance  for  the  two  populations.  The 
covariance  of  disturbance  events  at  patches  1  and  2  is  given  by  the  equation 


Cov(X1,X2)  =  E(XxX2)  -  E{Xx)  ■  E(X2).  (6.50) 

Only  one  of  the  four  possible  scenarios  ((1)  above)  results  in  a  11011-zero  value  for 
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covariance.  We  can  therefore  simplify  equation  6.50  in  terms  of  qi ,  p,  and  covariance 
of  disturbance  between  the  two  patches  (c)  to  for  our  purposes  to 

c  =  qi~p 2  (6.51) 

where  disturbance  tends  to  occur  at  both  patches  simultaneously  if  c  >  0,  and  tends 
to  occur  at  only  one  of  the  two  patches  at  a  time  if  c  <  0. 


6.6.2  Derivation  of  c  and  p  Inequalities 

If  c  =  q\  —  p2,  and  (6.45)  through  (6.51)  are  true,  then  the  following  inequalities  are 
also  true: 


Inequality  I 


Rearranging  (6.46), 
Qi=P~  92- 

Given  (6.51)  and  substituting, 
c  =  {p  ~  92)  -  P2- 
Rearranging, 
p-p2  -c  =  q2. 

Since  0  <  q<2  <  1, 

0  <  P  —  p2  —  c  <  1 
and  therefore 
c  <  p(l  -p). 
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Inequality  II 


Rearranging  (6.51), 
<?i  =  c  -  p2. 
Since  0  <  q\  <  1. 
0  <  c  —  p2  <  1 
and  therefore 
c  >  —  p2. 


Inequality  III 


Given  (6.45)  and  substituting  (6.49), 

9i  +  2q2  +  Q4  =  1. 

Rearranging  (6.46), 

92  =  P  ~  Qi  ■ 

Combining  these  gives 

9i  +  2  (p  —  91)  +  94  =  1 
and  2p  —  9!  +  94  =  1. 
Rearranging  6.51,  91  =  c  +  p2. 
Substitution  gives 
2 p  -  (c  +  p2)  +  9-1  =  1 

and  2/;  —  c  —  p1  —  1  =  94. 

Since  0  <  94  <  1, 

0  <  2p  —  c  —  p  2  —  1  <  1 . 
Rearranging, 
c  >  —p2  +  2p  —  1 
c  >  —  (p2  —  2p  +  1) 
c  >  —  (1  —  p)2. 
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Chapter  7 
Conclusions 

7.1  General  Summary 

The  goal  of  this  dissertation  was  to  explore  metapopulation  dynamics  and  population 
connectivity,  with  a  focus  on  the  softshell  clam,  Mya  arenaria.  I  first  worked  towards 
developing  a  method  for  using  elemental  signatures  retained  in  the  larval  shell  as 
a  tag  of  natal  habitat.  In  concert  with  this  method,  I  designed  and  implemented 
experiments  to  understand  environmental  and  biological  factors  that  may  influence 
elemental  signatures.  Next.  I  used  a  genetic  marker  to  study  regional  patterns  of 
diversity  over  the  current  distribution  of  M.  arenaria ,  including  populations  outside 
of  the  Northwest  Atlantic.  Finally,  1  used  theoretical  approaches  to  understand  the 
importance  of  migration  in  determining  patch  connectivity,  and  the  impacts  of  dis¬ 
turbance  on  met.apopulation  patch  dynamics. 

7.2  Chapter  Summaries  &  Future  Directions 

In  Chapter  2,  I  reported  on  experiments  designed  to  determine  whether  laser  abla¬ 
tion  ICP-MS  was  able  to  remove  and  measure  the  first  larval  shell  of  M.  arenaria 
reliably,  without  burning  through  to  the  juvenile  shell  material  underneath.  1  found 
that  ICP-MS,  in  its  current  configuration,  consistently  ablated  material  beneath  the 
first  larval  shell,  rendering  the  measurements  unreliable  for  determining  natal  habitat 
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based  on  elemental  composition.  Results  from  Chapter  2  have  consequences  for  in¬ 
vertebrate  ecologists  who  wish  to  use  elemental  chemistry  to  explore  questions  about 
population  connectivity.  Methods  that  have  proven  useful  for  other  study  systems, 
such  as  otoliths  in  fishes,  should  be  carefully  examined  and  tested  before  they  are 
transferred  to  a  new  system. 

Although  I  was  not  able  to  measure  connectivity  of  New  England  softshell  clam 
populations  using  their  larval  shell  elemental  composition,  there  may  be  other  sensors 
or  configurations  of  the  ICP-MS  that  would  be  able  to  accomplish  this  task.  For 
these  studies,  or  any  study  that  uses  shell  elemental  chemistry  as  an  environmental 
indicator,  it  is  important  to  understand  the  effects  of  factors  potentially  affecting 
shell  composition.  Such  factors  may  be  environmental  (e.g.  temperature,  salinity, 
seawater  elemental  concentrations)  or  biological  (e.g.  growth  rate,  age,  physiological 
condition). 

Chapter  3  investigated  the  effects  of  biological  factors  on  elemental  composition  in 
juvenile  M.  arenaria  shells,  specifically  age  and  growth  rate.  I  found  that  growth  rate 
was  negatively  correlated  to  some  elements  and  unresolved  in  others;  one  explanation 
of  the  negative  correlation  with  growth  rate  is  that  ability  to  discriminate  against 
non-Ca2+  ions  increases  with  the  animal’s  size.  I  found  that  age  significantly  affected 
incorporation  of  some  elements  into  shell.  This  was  evident  in  the  large  differences 
in  shell  chemistry  earlier  versus  later  in  the  juvenile  clam’s  life.  Both  of  these  results 
suggest  that  age  and  growth  rate  should  be  considered  when  designing  and  imple¬ 
menting  studies  where  shell  elemental  composition  is  used  as  a  proxy  for  environment, 
including  connectivity  studies. 

In  Chapter  4,  I  continued  with  my  investigation  of  the  factors  affecting  shell  el¬ 
emental  composition.  I  determined  the  effects  of  temperature  and  salinity  on  shell 
composition  of  larval  and  juvenile  M.  arenaria ,  and  compared  the  effects  of  these 
factors  between  the  two  stages.  I  found  that  both  of  these  environmental  variables 
influenced  shell  chemistry  for  some  elements  in  both  life  stages,  but  this  relationship 
was  different  for  the  two  stages.  The  conclusions  of  Chapter  4  are  therefore  similar  to 
those  of  Chapter  3:  environmental  factors  have  the  potential  to  affect  shell  composi- 
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tion,  and  therefore  must  be  taken  into  account  in  studies  that  use  biogenic  carbonate 
as  an  environmental  indicator.  It  is  important  to  note,  however,  that  geographic  vari¬ 
ability  in  temperature  and  salinity  may  serve  to  create  unique  elemental  signatures 
in  habitats  that  might  otherwise  produce  similar  shell  chemical  compositions. 

Chapters  3  and  4  are  useful  for  future  studies  that  use  elemental  composition  of 
shell  as  an  indicator  of  environment,  for  the  purposes  of  studying  population  connec¬ 
tivity,  climate,  pollution,  or  otherwise.  Experiments  like  those  reported  in  chapters 
3  and  4  are  essential  for  interpreting  the  results  of  fine-scale  chemical  analyses  that 
modern  technology,  such  as  ICP-MS,  has  made  possible. 

In  Chapter  5, 1  investigated  the  genetic  structure  of  M.  arenaria  populations  across 
the  current  distribution  of  the  species,  testing  for  patterns  of  regional  differentiation. 
I  found  that  populations  exhibited  extremely  low  genetic  variation,  with  a  single  hap- 
lotype  of  the  COI  gene  dominating  at  all  of  the  sites.  Populations  sampled  outside 
of  the  Northwest  Atlantic  exhibited  higher  diversities  despite  their  more  recent  in¬ 
troduction.  Chapter  5  has  implications  for  management  of  M.  arenaria.  One  of  the 
most  commonly  used  management  strategies  for  increasing  softshell  clam  populations 
in  New  England  is  by  seeding  less  productive  flats  with  juvenile  clams  that  usually 
originated  from  other  locations.  Our  results  suggest  that  this  management,  strategy 
is  not  likely  to  significantly  impact  the  regional  distribution  of  genetic  variation. 

I  developed  a  matrix  population  model  in  Chapter  6  to  explore  the  relationship 
between  migration  rates  and  metapopulation  growth  rate,  and  to  understand  the  ef¬ 
fects  of  disturbance  on  metapopulation  dynamics.  Although  my  model  was  applied 
to  M.  arenaria  for  demonstrative  purposes,  it  is  broadly  applicable  to  any  population 
that  experiences  stochastic  environmental  variability.  First.  I  used  the  model  to  de¬ 
termine  what  population  parameters  should  be  maximized  to  most  efficiently  increase 
metapopulation  growth  rate.  I  found  that  it  is  usually  best  to  focus  efforts  on  the 
patch  with  the  highest  individual  population  growth  rate.  It  is  also  important,  how¬ 
ever,  to  determine  the  spatial  distribution  of  individuals  in  the  metapopulation  (i.e. 
the  migration  rates).  If  migration  is  such  that  the  majority  of  individuals  are  located 
in  the  lower-quality  patch,  it  can  be  more  efficient  to  focus  efforts  on  protecting  that 
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patch  over  the  higher-quality  patch.  I  then  introduced  disturbance  and  reexamined 
the  results,  finding  that  they  held  true  when  the  environment  varied  stochastically. 
Finally,  I  expanded  the  model  to  include  two  stages  in  each  patch.  For  the  particular 
set  of  parameters  I  chose,  I  found  that  it  was  always  best  to  focus  on  the  higher 
quality  patch.  The  addition  of  the  second  stage  appeared  to  reduce  the  importance 
of  where  individuals  are  distributed  within  the  metapopulation. 

Metapopulations  are  a  useful  framework  under  which  to  consider  many  marine 
populations.  Coastal  benthic  invertebrates,  like  M.  arenaria,  are  especially  suited 
for  analysis  as  a  metapopulation:  they  occur  in  geographically  separated  habitats 
that  are  connected  via  dispersal  of  the  larval  phase  in  the  water  column.  Measuring 
larval  connectivity  (i.e.  migration  rates)  is  difficult  due  to  the  characteristics  of  the 
larval  phase,  but  this  fact  does  not  diminish  the  importance  of  migration  processes 
for  influencing  metapopulation  dynamics.  Therefore  it  is  important  to  have  both 
theoretical  studies  (e.g.  Chapter  6  of  this  dissertation)  that  quantify  the  importance 
of  migration,  as  well  as  empirical  studies  (e.g.  Chapters  2-5)  that  seek  to  measure 
actual  migration  rates. 
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